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1. Introduction. The numerical solution of partial differential equations by finite 
differences can be accomplished either by explicit or by implicit difference schemes, 
depending in great measure upon the class of problems to which a given equation belongs 
and also upon the subsidiary conditions which are imposed upon the solution. For 
instance, if a hyperbolic differential equation with Cauchy boundary conditions is under 
consideration, it is generally possible to find an explicit difference scheme according to 
which the solution at any given instant is determined explicitly in terms of known values 
at a few preceding time intervals. On the other hand, it is inherent in the nature of an 
equation of elliptic type, with Dirichlet data, that the solution at any interior point of 
the region throughout which the equation is to be satisfied depends upon the data at 
every point of the boundary of the region. Consequently, in this case, any explicit 
technique must involve all the values prescribed at the boundary. In general, even if 
such an explicit difference relation can be found, it is not in a form which lends itself to 
practical computation. It is usually necessary, therefore, to have recourse in this case 
to some implicit difference scheme, that is to say, one in which an undetermined value 
is defined always in terms of relations among other values, equally undetermined. The 
solution of such implicit schemes in most instances is extremely difficult, as it reduces 
ultimately to the solution of simultaneous equations involving a large number of un- 
knowns interrelated in a manner quite unsuitable for easy determination. If the system 
is linear, there are several classical methods by which it may be treated; most notably, 
by direct inversion of matrices of large size [1], by iterative procedures [2] and by relax- 
ation [3]. It is the purpose of this paper to introduce a finite difference technique especially 
adapted to the solution of boundary value problems. This method depends upon the 
selection of an implicit difference scheme involving relations among only a limited 
number of the unknowns at a time, in a form that can be conveniently solved. It is 
especially well suited for use with digital computing machinery since the demand for 
memory capacity is small and the time required for computation is short compared with 
other existing methods, even when meshes of large size are involved. 

In this preliminary report, the conditions for stability of the process are investigated. 
Criteria are derived for stars appropriate to various classes of partial differential equa- 
tions. Applications of the theory to problems of both hyperbolic and elliptic type are 
discussed and the advantages of the method are illustrated. 
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Fig. 1. Typical meshes for boundary contraction 


2. Boundary contraction. Consider a region ® (Fig. 1) bounded by a piecewise 
smooth Jordan curve S, and let ® be referred to an (r, @) coordinate system topologically 
equivalent to the polar system with the contour S, taken to be the unit circle. A sequence 
of concentric contours S, , S, , --- of decreasing radius, which will be referred to as circles 
in the following, may be defined interior to S, and a mesh fitted onto ® by considering, 
as nodal points of the mesh, the points of intersection of the circles of this sequence 
with a finite number (N + 1) of radii from the origin. Suppose that the prescribed data 
are sufficient to determine initially the values of the solution at the nodal points on, 
say, (j + 1) of the outer circles S, , S, , --- S, . For instance, if Dirichlet data are given, 
S, is prescribed initially; if Cauchy type data are prescribed involving both the values of 
the solution as well as its normal derivative at the outer edge, then S, and S, are given 
initially. The contraction method for solving boundary value problems is a step-by-step 
process that is initiated at the outer boundary, where the originally prescribed data are 
given. At each stage, the solution is sought only at the nodal points corresponding to the 
outermost circle where the solution remains unknown. The newly computed data on this 
circle then replace those used in determining it and a new boundary value problem, 
similar to the original one, arises on a somewhat smaller region interior to the previous 
one. For example, if the data are known on S, , S; --- S; , this information is used to 
compute only the unknown data on S;., ; then, the known information on S, , S.,---, 
S;.: is used to determine S;,.. ; and so forth. The new problem is treated in a similar 
fashion and the solution moves inwards progressively with the boundary being con- 
tracted onto the center. It should be observed that the data, as specified at the outset, 
are soon forgotten but that the information is propagated to the center in modified form 
as the new conditions on the inner boundaries. Thus, when using this method, the tre- 
mendous advantage is obtained that at any given stage it is necessary to work with 
only a small proportion of the total number of nodal points of the mesh. 

3. Stability and propagation of error. The questions of stability and propagation of 
error are critical ones for this numerical method. This may be understood most easily 
in terms of the notations that will be used subsequently. Let the angular arguments of 
the (NV + 1) radii, to which reference was previously made, be denoted as @,, (n = 0, 
1, 2, --- , N) and u(k, n) = u(r, , 6,) be the value of the solution at the nodal point 
defined as the intersection of the radial line determined by 86, with the circle S, , (k = 0, 
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1, 2, ---) which has radius r, . Also, let v, be the column vector with the components: 
u(ke, 0), u(k, 1), «++ , u(k, N) which are the values of the solution on the circle S, . The 
contraction process consists in the repeated application of a transformation A, that 
relates the (j + 1) vectors Vy4j41 , Ve+j » °°" » Ve+s to the (j + 1) vectors Vy.; , Ve+j-1, 

, v, . In the particular case of linear differential equations to be considered here, 
the transformation is expressible in matrix form. Using block representation, this can 


be written in the form: 


Vivi | | Vi 
Tare | Vis 
= A, (1) 
Visjas Vis; 
Thus 
{VierVere °°° Verses = [] A. {vov, 22° Vit. (2) 
+=0 
or typographical reasons, here and subsequently, it is convenient to write column 


vectors row-wise, using braces, as: 


u(O, 1) | {v,| 
u(O, 2) | Vv; 

= {u(0, 1)u(O, 2) --- u(O, N)}; - | = {vov, «+ vj}. 
u(Q, N)] ly, 


To illustrate the notation with an example, suppose that an explicit star relation is 
given by: u(2 + k,n) = 1/2[ul + k,n — 1) + u(l + k,n + 1)] + 1/4 [u(k,n — 2) + 
Qu(k, n) + u(k, n + 2)] where the indices corresponding to the angular argument are 
assumed to be reduced modulo N + 1 and k = 0, 1, 2, --- . Then: v..2 = Pv,.,; + Qv, 
where P and Q are the (V + 1) X (N + 1) circulant matrices 








01 0 0 1| (210 -+ O 1 
pai|! 0 1 0 0) 9 =:/! 2 1 0 0. 
= . 

l1 0 0 1 0) lt 0 0 1 2) 


Using block representation: 
Vivi! _ ° I| iv, 
\P Q Views 





Vi+2 


where 0 and I are respectively the (V + 1) X (N + 1) zero and identity matrices. In 
this simple case, the matrix 


0 I 
P Q 


A= 
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is independent of k, so that 


Pa 
j ett) -_ Dead bd 


Vi Vi 


As k increases and the solution moves inwards, an error made at some point in the 
calculation may be propagated forwards and may give rise to a first type of instability. 
This may be examined by assuming that for some integer m, an error ¢,,,; is introduced 
in v,,., , and then determining the resulting error at a later stage. From (1) this is given 


DY fn42+)+1 , Where: 


m+k 
eo cacs °** Cesar euskesers = I] A;{0 --- O¢,,.;}. (3) 
‘=m 
The magnitude of the resulting error can be measured by comparing |e,,.,.;.: | With 
| ¢,,.; | Where | e | is the maximum modulus of the components of e. If v(r) is the vector 


consisting of the values of u on the mesh circle of radius r, and if an error occurs in v(r’) 
for some r’ > r, then the resulting error in v(r) will depend upon the choice of mesh used 
in going from 7’ to r, since, in general, the matrix A; is determined by the mesh. For 

ach mesh, let m + j be such that v,,.; = v(r’) and k such that v,,...;+1 = v(r). Assuming 


that the solution is uniformly bounded over @&, for stability of the contraction process, 
it is required that « ;+, be uniformly bounded as the mesh decreases, independently 
of r and r’. This in turn, requires that the matrix product ([]‘.. A.) should remain 
bounded for all k as the mesh decreases. Thus it is possible, in further consideration of 
the first type of instability, to pee discussions to this matrix product. 

It should be noted en at Eq. (1) is an explicit vector relationship to determine v,.;+1 ; 
however, the components of i. vector need not be explicitly determined ay they are, 
indeed, usually ad implicitly one with another. The solution of these implicit 

elations gives rise to a second type of instability which is described best, perhaps, by 
the aid of a simple example. In particular, for the case k = 0, j = 0, consider the approxi- 
mating star relating the components of v, to those of vp : 


Cc. Uln +2) +e, wl,nt y+ e6,u(1,n) +e, uO, n + 2) (4) 
+¢°u(0,n +1) +c, u(0,n) = 0 (n = 0,1,2,--* , AN 
where the indices corresponding to the angular argument are assumed to be reduced 


modulo (N + 1). If c;’? ¥ 0, this may be solved as: 


u(1,n + 2) = _(- : ) > ¢ ul,n +8) — (. ) Dd ofulO, n + 8), (5) 
‘d — cs — 


1, 2, --» , N, ul, n + 2) is computed 


which implies that, for each n = QO, 
from u(1, m + 1), u(1, mn) and the nodal values on S, . The system (5) can be expressed 


as a matrix equation: 
wi ; , 1 , e 
u(l,n + 1) ul,n + 2)} = Giu(l,n) ul,n+1)} - a H,R’v, , (6) 


where: 
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, . (7) 


(1) (1) (1) (1) 
—Co /c2 =Z /C2 


G = 


H,, is the 2 & (N + 1) matrix: 
0 0 0 O --+- O 


H, = | (8) 
oe ag” og” 0 +++ 0 
and R is the (V + 1) & (N + 1) matrix: 
0100. 0 0 
0010. 0 0 
oOo 81+ @ @ 
(9) 


oOo 0 @ ~- i © 
1000. 0 0 








The system is seen to be in a very convenient form for numerical solution since Vp is 
already known, so that once u(1, 0) and u(1, 1) are determined, then u(1, 2) can be 
computed explicitly from (5) with n = 0; similarly, with w(1, 1) and u(1, 2) not known, 
u(1, 3) follows with n = 1; and so on. To determine u(1, 0) and u(1, 1) it is merely 
necessary to start the process with a good initial guess for these quantities and perform 
the indicated calculations by making a complete circuit around S, , coming back to 
, 0) and u(1, 1). In virtue of the continuity of u in @, the results obtained at the end 
of the circuit will agree with those at the beginning, whenever (4) is a linearly independent 
system of equations, if and only if they are the correct solutions. Thus, the implicit 
scheme (5) really reduces to two simultaneous equations in two unknowns. Analytically, 
the system (6) may be solved in the same way by carrying out the indicated substitutions 
and equating the end result to the initial value. Then: 


fu(l, N + 1) ul, N + 2)} = {u(1, 0) u(1, 1)} 
. 1 NV : (10) 
GY **{u(1, 0) u(1, I} — (-1;) >> G’"'HR'y, . 
It is clear from (10) that the numerical calculations will become unstable as they proceed 
around S, , if the iterated matrix G is unstable. Thus, a second stability condition, in 
addition to the one previously indicated, must be satisfied if a boundary value problem 
is to be solved satisfactorily by the contraction method. 

4. Instability of the first type. In this section criteria wil] be developed for stability 
of the first kind discussed in Sec. 3 above. At the outset, a restriction is imposed upon 
the type of star to be considered and it will be assumed that this is limited to an approxi- 
mating function of the general form: 


N 
Ye uk + jti,ntsy + DcPuk+jn+s+-:- (11) 


e=U0 
Vv 


> cfu(k,n + 8) = 0, (n = 0,1,2,--- ,N;k = 0,1, 2, -->). 


— 
Ppa 


It should be noted that the weights c‘”’ in this svstem of equations are independent of 


n as well as k, which implies that the basic form of this star is unchanged as contraction 
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moves inwards and also as the computation moves around S,,;,; . A star of this type is 
appropriate, for instance, to homogeneous linear partial differential equations which 
have constant coefficients when referred to the chosen mesh system. Certain generali- 
zations upon this will be indicated later in Sec. 8. 

The stability criterion for the system (11) corresponds in a simple way to the system 
itself. To illustrate this, the final results obtained in this section will be stated here. 
By taking k = n = 0, a set of (N + 1) algebraic equations can be derived from (11) by 
replacing the «(0, s), w(1, s), --- , u(j + 1, s) with w (where w is any of the (N + 1) 
roots of w** = 1) and multiplying each sum starting from the left by \’*', \’, --+ , A, 
1 to obtain: 

N N V N 
NS co HAM Velo +e +A Vela’ + Deeiw’ = 0. — (12) 
s=0 2=0 s=0 s=0 
The essential condition for the stability of (11) as contraction moves inwards relates 
to the roots of (12). If the circulant matrix C‘’*"’ defined by (14) is non-singular, necessary 
and sufficient conditions for the stability of the system (11) are that for each value of w, the 
(j + 1) roots of (12) should be distinct and have moduli less than or equal to unity. 

It may be remarked that in practice this stability criterion reduces to an exceedingly 
simple condition since a great number of the weight factors are usually zero. For instance, 
the stability equation associated with (4) is: 


Mes?a? + efPw + of?) + Cw? +528 +5) = 0. 
For the derivation of this result, Eqs. (11) can be combined into a matrix equation: 
ee vss FO ay Hoe + OO + Cv, = 0, (k=90,1,2,---), (13) 


where C‘” is the circulant matrix: 


(i) i i) i) 
Co C; C2 * Cy 
(¢) (#) i) (s) 
\Cn Co C; ° Cy-1 N 
i) (i) (i) i) i) ( 
C = |Cy-1 Cy Co * Cy-2| = er. (14) 
. . | o=s 
| 
( ( +) a) | 
le," C,"” C3 * €6 


If Eqs. (11) are sufficient to determine u(k + j + 1, 0), u(k + j + 1, D), 
-++,u(k + 37+ 1, N), the system must represent a linearly independent set of equations 
in these unknowns; this implies that C‘’*’’ must possess an inverse (C“*")~" = D and: 


Vi+j41 = — D(C“ v,,; +C St Re + 5 pes + C’y,). (15) 
It is convenient to express (15) in a form similar to (1): 
{Vi+iVe+2 °° Veaser} = AlViViss -* Vest (16) 
with 
[ o I 0 - 0 
.*s 0 I 0 
A= | , ‘ ; ; . : (17) 


| 0 0 0 : I 
(-pc® -—pc” -pc” - -—DC” 


where I is the (V + 1) X (N + 1) identity submatrix. It should be observed that by 
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virtue of the assumption that the weights c‘"’ do not vary with r, the matrix A does not 
change as contraction proceeds. Equation (2), therefore, in this case becomes simply: 


fviaiVaee °° ° Vassar} = A’ {vov, --- V;}. (18) 


It is familiarly known [4] that the stability of an iterated matrix depends upon its 
eigenvalues and eigenvectors and that if the matrix has a full complement of eigenvectors 
then a necessary and sufficient condition for its powers to remain bounded is that the 
moduli of its eigenvalues are less than or equal to unity. The eigenvalues and eigenvectors 
of A can be found, by virtue of the fact that the matrices C“” are circulant, or rather, 
by virtue of the fact that the weights cS” do not vary with 6. Let \ be an eigenvalue of 














A and {x,, x, , :-+ x,} an eigenvector, where each x, has the same dimension as v, . The 
relation 
0 I . 0 ] (x, x. | 
0 0 ° 0 xX; xX, 
0 0 . x X2 ; 
re |e ans (19) 
0 0 . I X;-1 X;-. 
—pce® -pce” - —DC”) lx; | (x; 
gives X;., = Ax, (¢ = 0,1, --- , 7 — 1) so that the eigenvectors of A are of the form: 
{xX AXo A°*Xy +++ A’Ko}. (20) 


The vector x, satisfies 
afer’ 22" 4 ee 4. 3s ee 4 Cr, aw GS, (21) 


which can be written as 


N NV N N 
(, TS cl tPR' + SeclPR' + --- $A DocMPR’ + > o!°R’)x, =0, (22) 
* ) e=0 s=0 s=0 
where R is defined by Eq. (9). The eigenvectors of R are {1w w +++ w'}, where w 
ranges over the (NV + 1) roots of w* = 1. If x) = {l ww’ --- w’}, Eq. (22) gives: 


N N N N 
Ne De Pat HM Dice tee FA Dice + Dieiw = 0. — (28) 

s=0 2=0 s=0 
For each value of w, (j + 1) values of \ are determined and if these are all distinct, (j + 1) 
distinct eigenvectors of A are found from (20). As w ranges over the (N + 1) roots of 
unity the full complement of eigenvectors of A are obtained and A is stable, provided 
that all the roots of Eq. (23) have moduli which are less than or equal to unity. This is 

the result stated at the beginning of the section. 

It is interesting to comment upon the type of instability which occurs when the 
above conditions are relaxed. This is most easily discussed by consideration of the 
classical canonical form B to which A can be reduced by a similarity transformation 
[5]. The matrix B is composed of a number of primitive blocks (8;) of the type 


(a. 10. 0 

Oo & i « ® 

ict. « et %% % (24) 
lo 00-1 

eoee8- 
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Since A is similar to B, the powers of A are similar to the powers of B and the behavior 
of the latter is as the behavior of the powers of its primitive blocks (8;). In case any of 
the roots \ of the system of Eqs. (23) has modulus greater than unity, then B* increases 
exponentially like \* so that B, and consequently A, is unstable. In case all the roots of 
(23) have moduli less than or equal to unity but one or more has multiplicity two, then, 
at worst, B* increases like 


i, ] = O(k). (25) 
0X 


More generally, if one of the roots has multiplicity m, then B* = O(k"~'). In such in- 
stances, the matrix is subject to a controlled instablity, which may be suitable for 
many problems when the number of contracting steps inwards is not large. 

5. Instability of the second type. In the example of Sec. 3, illustrating the second 
type of instability, Eq. (4) was solved for u(1, n + 2). Provided cj’ # 0, the same 
star could be used to compute u(1, n) from u(1, m + 2) and u(1, n + 1); or alternatively, 
if c’” 0, to compute u(1, m + 1) in terms of u(1, n + 2) and w(1, n). In each case, 
however, a different stability condition is implied. The generalization to the more complex 
star (11) is quite simple but involves some notational inconveniences. Suppose that 
cS'*? = 0 and for each n = 0, 1, 2, --- , N the value u(k + j + 1, + p) in (11) is 


computed from 


, 


i+1), - . 
cc uk+j7+1,n+58) 
sa (26) 


— 
— (1/e8'*?) SS Soc ulk + i,n +), 
t=0 s ) 


uk+j+i,n+p) = —(1/c{'*”) 


N 


where the prime on the sum means that s + p. 

It is implied by (26) that an initial guess for v,,;,, is made, u(k + j + 1, n + p) 
is evaluated from Eq. (26) and then this newly computed value replaces the previous 
estimate. The same procedure is applied to each mesh point of S,.;,, in turn and at the 
end of the first circuit the new values are compared with the old. If they do not agree, a 
second circuit is made and the process is continued until the solution of the system is 
obtained. With the indices corresponding to the angular argument reduced modulo 
(N + 1), Eq. (26) can be written in matrix form as: 


(uk +j+i,n+2+>) uk+jt+in+1+>p)| 
luk +j7+1,n+3 +p)! uUk+j7+1,n+2+p) 
uk + i+ in+4+p)| = Glue +it+in+3 +p) (27) 
luk + 5+ 1, n + p) uk+j+iln—1+>p) 


i 
— (1/c)’*”) > HR’v,..; , 


s=0 


where 
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> 








| 0 1 . 0 
0 0 . 0 
0 . 0 
G=a| | - + (28) 
—ciit” Jes (7+1) —cii3” i , —cfiy? ‘oes 


and H; is the N X (N + 1) matrix: 


0 0 0 0 
0 0 0 - O 
H; as ° ° ° ° ° (29) 
0 0 0 0 
| (4) (#) (#) (#) 


Co G C2 * Cy 


From Eq. (27), at the end of a circuit: 


+ j+i,n+1+>p) uk+j+i1,n+1+>p) 
G41 n-+2+ p) — Gg sili Nici iad lid (30) 
lu(k + 7+ La-tee uk-+j+i1,n—1+p) 
nate (1/cS*”) > > G*"'H.R**'v,.; : 
i=0 t=O 


It is clear from Eq. (30) that in making the circuit around S,,;,, any error which may be 
introduced will be propagated in a manner which depends upon the properties of the 
iterated matrix G. To guarantee stability of G, it is essential that its non-zero eigen- 
values should be distinct and have moduli less than or equal to unity. But G is seen to 
be in its Jordan canonical form [5] so that its characteristic equation is simply: 


i 
; 1 N . i+1) N-1 3+1) N-2 (j+1) (7+1) 
C, MepCp-1 ob + Cp-2 b $+ ose + Coes K + Cp+1 = 0. (31) 


Thus, the second stability criterion is that the non-zero roots of (31) should be distinct and 
should have moduli less than or equal to unity. 

If the solution of Eq. (27) is to be found iteratively by using the results after one 
circuit as the initial guess for the next and if this process is to converge, conditions 
slightly stronger than those for stability must be satisfied. Jt is essential for convergence 
that the moduli of the roots of (31) should be less than unity. 

It will be recalled that the first stability conditions required that the circulant matrix 
C‘’*” should have an inverse. This is equivalent to the condition that it has no zero 
eigenvalues. If C‘’*” is written: 


N 
cut a } ee (32) 
s2=0 


then its eigenvalues are given by: 


N 
y die, (33) 


a=0 
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where w*** = 1. This is equivalent to: 
N 


N N 
+ ey = bE nde x p—s) ‘ate Pitas awe, (34) 
s=0 s=0 s=0 
where the subscript indices of the weight factors are reduced modulo N + 1. From 
this it follows that a sufficient condition for C’*" to be regular is that the equation 


N ° 
> ci yy? = 0 (35) 


should have no roots on the unit circle. Identifying this with (31) it is seen that the 
criterion for convergence is sufficient to ensure the regularity of C“’*”’. 


6. Example of an equation of hyperbolic type. It is required to find the solution: 


u = u(r, 6), of the hyperbolic equation r° (d°u/dr*) — (d°u/d6") = O in the region 
0 <r Sl, satisfying u(1, 6) = f,(@), du(l, 6)/dr = f.(0) on the boundary 


and | u(r, @) | S M in the interior of this region, where M is some constant. As mesh, let 
S, be the unit circle, S, the circle of radius r, and choose r,,,/r, = p where p is a constant 
0 < p < 1;also, let u(k, n) = u(r,, NAO), where A@6 = 2r/(N + 1) for some integer N. 


The first boundary conditions are satisfied by taking u(0, n) = f,(n A@). The de- 
rivative appearing in the second boundary condition may be replaced by its divided 
difference approximation: du(1, 6)/dr  [u(O, n) — u(1, n)]/(1 — p) so that: u(l, n) = 


u(0, n) — (1 — p)f2(nAé). 
It is possible to select a suitable star for this problem that involves eleven nodal 
points so that only ten multiplications are necessary to compute one of the nodal values 


from the rest. A typical stencil is given in Fig. 2. 


4 


k+l, m2 Sok, net 








k+1, n=2 


Fig. 2. Stencil for hyperbolic equation 


In this case, the solution on S,,. is determined always in terms of its known values 
on S,.;, and S, so that j + 1 = 2 in (11). 

An approximation for the second derivatives appearing in the differential equation 
may be defined by taking linear combinations of the fundamental difference approxi- 
mations at different points in the stencil. Thus for any selection of the real numbers a 
b, c such that a + b + ¢ = 1, a difference approximation may be defined 


ui 1 pag —e sisi 
ag ~~ ae fafu(k + 1,n + 2) — Quk +1,n+1) + uk +1,n)] 


+ bluk +1,n+ 1) — 2uk+1,n) +uk+1,n — 1] (36) 
+ clu(k + 1,n) — 2uk+1,n —1) +u(k+1,n — 2)]}. 
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In a similar way, using divided differences: 


, OU = 
Ss i ee : “ 
* or pi 1 -— p) | | + p) 


+ pu(k,n + 1)] + efu(k + 2,n) — (1 + p)u(k + 1,7) + pu(k, n)] (37) 
+ flu(k + 2,n — 1) — (1 + pu(k + 1, — 1) + pulk,n — 1)]}, 


{dlu(k + 2,n + 1) — (1 + p)u(k + 1, + 1) 


where d + e + f = 1. The stability criteria may be satisfied for a wide range of values 
of the parameters a, : c, d, e, f and this fact may be used in actual practice to assist in 
the selection of an approximating star which is the most accurate possible, compatible 
with these condit : ms. However, we are not here concerned with the problem of accuracy 
and we therefore assign the parameters in such a way that the resulting algebra is some- 





what simplified. Choosing a = 1,b = — 1/4,¢ = 1/4;d = l,e = — 1/4, f = 1/4, the 
following difference approximation to r* 0°u/dr? — d°u/d6? = O results: 
9 
= ee fu(k +2,n+ 1) — 4u(k + 2,n) + 4u(k + 2,n — 1)] 
P F aa 
SS +2) + (~ - — 2 + 1,2 + 1) 
= Ae p(l— ps 
; - 
—; = k 38 
ai fia * 4A0 mr oe i. , n) (38) 
+ (_3 se dpi ia s)u(k +1i,n-—1) - “ne u(k + 1, — 2) 
GAG p(l — p)/ , 4A 
9 
+ G— i+ - ~ [uth ,n + 1) — Rutk,n) + fu(k,n — 1)] = 
From Eq. (12) the first stability equation is 
4 4 ‘ a as il — ew 
(1 —p(l+op Fez + p(l — p) 4A" - 
‘.. — 7os)o" sl ot 
+(z LAG 2p(1 — =) ¥ (si — p)” 4A + rr ae wh (39) 
4 Aw — i tio) _ Ww —- 7+ we) {n° 
"ai-pits+p Ali-pi¢+a 
p(l — p) Ul + p) (. 2 - 3) . " 
2 ae Tne ee ae” Pe 
where w = exp {2nzi/(N + 1)}, (n = 0,1, --- , N). It is advantageous geometrically 
to take A@ = (1 — p), so that Eq. (39) reduces to: 
 — (1 + al — p) + pcos (<% Z =)h + p= (40) 


A quadratic equation: x* + px + q = 0 with real coefficients has roots with moduli less 
than or equal to unity in case 1 = g 2 | p| — 1. Application of this result easily demon- 
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strates that the roots of Eq. (40) have moduli less than or equal to unity. The equation 





> — (1 + p)[(1 — p) + p cosg]A + p = 0, (41) 
has equal roots only if 
p+ 2(p)* — 1 
cos @ = ; (42) 
¢ p(1 + p) 


and Eq. (40) will have equal roots only if ¢ = 2nx7/(N + 1) forsomen = 0,1, 2, --- , N. 
Since Eq. (42) can be satisfied, if at all, only for isolated values of n, it is not difficult 
to select values of N so that the roots of Eq. (40) are distinct. 

The second stability equation, from (38), is seen to be 


up —-ipt+}=0, (43) 


and it is readily verified that the roots of this equation are distinct with moduli less 
than unity. 

With the boundary conditions taken to be f, = cos 0, f. = 1/2 cos 6, the analytic 
solution: u = r’”’ cos [(3)'”” (log r)/2] cos @ for the differential equation can be found. A 
comparison between the computed and the actual values indicates that with (V + 1) 
as great as 1000, stability was indicated for as many as 250 steps inwards so that a mesh 
of 250,000 nodal points was involved. 

7. Example of an equation of elliptic type. The most important application of the 
boundary contraction method is to partial differential equations of elliptic type. No 
essential difficulty exists as a consequence of the well-known distinction between elliptic 
and other types of partial differential equations. In practice, however, it is not easy to 
find difference schemes for this class of equations that fulfill the stability requirements. 
In another paper [6], the authors have presented a detailed analysis of Laplace’s equation 
in the circle. It was found that for this equation a modified version of the method pre- 
sented here was preferable. However, the stability criteria related back to a simple form 
of those which have been given here. Solutions have been computed over a mesh of 2,500 
points and accuracy to six significant figures has been obtained for a variety of boundary 
conditions. Stability was verified up to 200,000 nodal points but similar accuracy cannot 
be claimed for meshes of this size. 


8. More general forms of stars. 


By arguments analogous to those used in the fore- 
going paragraphs, the contraction method may be shown to be valid equally well for 


stars of the form 


I : | 


N 
De Puk + Gti nta + VeMuk +jn+s + 
+ Dd cf u(k, n-+s) = F(r, 0 


which correspond to a non-homogeneous linear partial differential equation. Equation 
(44) can be expressed as: 
( 15 


} + 2) 


SVeaiViso °° * Veosei} AlviVie: °° Vios} +B, 
where B, is a matrix independent of the v, . It follows readily that the stability properties 
of Eq. (45) depend only on those of the matrix A, which has been fully discussed. 

A further generalization is possible to a system in which the weights c}'’ may vary 














1960] NUMERICAL SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS 13 


with the radius, as the contraction moves inwards. This implies from Eq. (1) that the 
matrices A, depend upon r. It is found that sufficient conditions for stability of the 


first type for a star such as: 


‘ 
Ye rosdlk +i ti,nts + Del Om. )uk +j,n+s +>: 
| _ (46) 


N 
+ > cf (r,u(k,n + 8) = F(r, 8) 
s=0 
are: (a) that the mesh be taken sufficiently fine, and (b) for each w, such that 
= 1, the roots of the equations 


— 


N N 
YEP Osea” HM DY Py jo’ tees + Yk r,)w* = 0 (47) 

e=0 s=0 
should be distinct and should have moduli less than unity, uniformly in k. This result 
follows from the proposition [7] that if a matrix A = A(r) is continuously dependent 
upon a parameter r which takes on a sequence of values: {r,; | 7 = 0, 1, 2, ---} in a closed 
bounded region, then the matrix product [], A(r;) is bounded, provided A(r) has a full 
completement of eigenvectors, the moduli of the eigenvalues of A(r) are uniformly less 
than unity and 6 = max, [|r;,, — 7,|] is sufficiently small. The second stability criterion 
for Eq. (46) becomes simply that the non-zero roots of 


Yc G sede = 0 (48) 


n=0 


should be distinct and have moduli less than unity. 

The generalization to stars in which the weights c}'’ are dependent upon the angular 
argument has not yet been made. 

9. Further investigations. [Further investigations of the contraction method are 
being made by the authors. The method is being generalized to include other classes of 
partial differential equations and to stars which depend on the angular as well as the 
radial arguments. The possible applications of the method to irregularly shaped bound- 
aries, to multiply connected regions, to mixed boundary value problems and to non- 
linear equations will also be considered. 
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BOOK REVIEWS 


Dynamic programming. By Richard Bellman. Princeton University Press, New Jersey, 
1957. xxv + 340 pp. $6.75. 


From a formal point of view, Dynamic Programming is the study of functional equations or of 
recursively defined systems of functional equations which involve a maximum operator. Its principal 
application is to the analysis of decisions over time or ‘‘multi-stage decision processes’’. Practically 
a single-handed creation of Richard Bellman, this novel technique has had a number of expositions in 
journal articles and shorter RAND publications. But the present book ‘Dynamic Programming’’, now 
in its second printing, stands out as the first systematic treatment at book length of this important 
and fruitful subject. It is an original and authoritative, attractive, sometimes brilliant, occasionally 
drawn out and repetitive, but throughout fascinating presentation of this lively and versatile mathe- 
matical tool. 

The general plan of the book is as follows. The basic ideas are introduced in terms of a specific 
problem concerning the repeated allocation of a single resource, an example that may well represent 
the simplest non-trivial case of a Dynamic Program. Following a second more complicated example 
involving probability (see below) the general structure of Dynamic Programming is exhibited, existence 
and uniqueness theorems are developed, the sensitivity of solutions and such properties of the solutions 
as convexity and continuity are discussed. An important tool in these investigations is the method of 
successive approximation through iteration. It is shown in particular that a suitable choice of the initial 
approximation will always result in monotone convergence. The remaining part of the book studies 
in detail various cases of a more advanced nature involving respectively, integral operators, systems 
of differential equations, partial differential equations, Min Max operators and Markov chains. 

Apart from the opening problem which has little intrinsic interest, the applications developed at 
length are essentially of three kinds. The first of these is illustrated by the so-called gold mining problem,— 
obviously a military mission problem in disguise-where the object is the repeated assignment of a 
given piece of equipment to alternative uses subject to constant (but different) probabilities of ruin 
and to diminishing returns. This problem is studied in both a discrete and a continuous version (chs. 
II and VIII) for two and for more alternatives, and its very simple and reasonable solution is given 
in explicit form. 

The second kind of problem is illustrated by a ‘‘bottleneck problem” (chs. VI, VII) known in the 
economic literature as a “dynamic input-output model’’, which involves the allocation of products as 
inputs to both current production and the expansion of productive capacity. As in the more conventional 
linear programming treatment of this model, advantage is taken of the existence of dual variables or 
“efficiency prices’”’ to obtain a sufficient characterization of the solution in terms of profitability con- 
ditions. However, the treatment is bogged down by a somewhat unfortunate choice of notation and 
never gets to the point of discovering in its general form the economically important principle of balanced 
growth. It is even doubtful whether the trial and error approach given here can be dignified by the name 
of Dynamic Programming. 

Both the goldmining and the interindustry example concern closed systems. Of much greater po- 
tential usefulness in economic applications seem to be cases of open systems, exemplified by the third 
main application in the book, the so-called inventory problem (ch. V). This is a classical piece of opera- 
tions research concerning the reordering of stocks by a firm facing a stochastic demand. Only the simplest 
case (linear ordering cost) is presented in full but several extensions are outlined. Recent work on the 
more interesting case which involves fixed ordering costs and especially some results by Scarf proving 
the optimality of the so-called s, S policy under very weak assumptions, have now rendered the state- 
ment obsolete (p. 153) that ‘‘at the present time practically no solutions of the corresponding functional 
equations exist and very little seems to be known concerning the character of the optimal policies 
arising from processes of this more realistic type.”’ 

As a demonstration of the theoretical power of Dynamic Programming, the book includes three 
interesting applications to mathematics proper: to the calculus of variations with special reference to 
cases involving inequalities as constraints (ch. IX); to sequential games, in particular to so-called 
games of survival (ch. X); and to the study of Markov chains (ch. XI). In ch. IX a set of partial differ- 
ential equations is obtained representing the more conventional Euler-Lagrange equations in an inte- 


(Continued on p. 80) 
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WAVE REFRACTION AT AN INTERFACE* 


BY 
C. M. ABLOW 
Stanford Research Institute, Menlo Park, Calif. 


Abstract. A plane wave in one of two perfect gases moves toward the parallel plane 
interface between the gases. The wave is either continuous or headed by a shock front 
weak enough that entropy changes may be neglected. Using Riemann’s solution of the 
appropriate hyperbolic partial differential equation, four equations are derived giving 
the details of the reflected and refracted wave motions. The equations are of first order 
integro-differential or implicit functional form depending on the boundary conditions 
and must be solved simultaneously for four functions of a single independent variable. 
The equations are suitable for numerical step-by-step solution. 

I. Introduction. The refraction of a non-uniform wave at a fluid interface is obtained 
by solving a second order partial differential equation under conditions along boundaries 
whose locations are not known in advance. Heretofore a finite difference numerical 
scheme covering a two-dimensional net has been used to establish the location of the 
boundaries and solve the boundary value problem simultaneously. (See, for example, 
Courant and Friedrichs [1], Sec. 83.) This report describes an analytic method reducing 
the refraction problem to the solution of one-dimensional integro-differential equations. 
The reduction in dimension is the principal advantage of the analysis. The analytic 
form may also provide a convenient starting point for existence or uniqueness proofs. 
The method is an application of Riemann’s solution of the partial differential equation 
involved (Courant and Hilbert [2], vol. II, Chap. 5, Sec. 4). 

The mathematical formulation of a number of interface problems, including those 
appearing here, is given by Drummond [3]. The application of this fluid theory to the 
refraction of shocks in metal is discussed in Duvall and Zwolinski [4], Walsh and Christian 
[7] and Walsh, Shreffler and Willig [8]. An analysis similar to the present one is given by 
Heinz [9] for the reflection of a shock wave in air at a rigid wall. 

II. Equations of motion. Each of the two fluids is supposed an inviscid, polytropic 
gas. It is further assumed that any shock discontinuities occurring are so weak that the 
entropy changes produced may be neglected. The initially homentropic fluids will then 
remain so. 

With these assumptions the equations of one-dimensional flow of either fluid may be 
written in characteristic form (Ref. [1], Sec. 37). There are two families of characteristic 
curves which are referred to as the C+ and C— characteristics in the (z, ¢) or physical 
plane and as the [+ and I'— characteristics in the (l, u) or hodograph plane. Along 
each characteristic of one family r is a constant, C-+:2, = (u+c)t,,andT+:l+u = 
2r. Along each characteristic of the other family s is a constant, C—:z, = (u — e)t, , 
and [—:l1 — u = 2s. Here 





2 @ 2 
= v - . ao = 
p=ap’-B8, =P, tae, 
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literal subscripts mean partial differentiation, and the letters have their usual meanings: 


distance from a fixed plane parallel to the interface, 


b = 

t = time, 

u = flow speed in z direction, 

c = local speed of sound, 

r, 8 = Riemann invariants, 

l = convenient thermodynamic function, 
? = pressure, 

p = density, and 

a, 8, y = given constants of the fluid (y > 1). 


The characteristic equations in the physical plane may be written 


y¥+1 ¥—3 
(pes 38) 


-_ 


, y¥—3 r+1,) 
On 3 ge | eee ce see 
° ( 5) r+ >) & }t, 


II 


C+ :2, 


or, eliminating z, 





t+. _ o. ne 
(A) roe. eS 8D 


If ¢ is found as a function of r and s satisfying equation (A), then the characteristic 
relations C+ and C— determine x as a function of r and s. Inverting these relations, r 
and s are then known at each point (z, ¢) in the physical plane or, using the [+ and 
I'— relations, u and J are known there. Since / is a univalent function or p or p, knowing 
any one of 1, p, p, or c at a point determines the other three. 

The solution of equation (A) under quite general boundary conditions has been given 
by Riemann. For that solution one uses Riemann function R = R(r, s; &, ) satisfying 
the following conditions: 

(a) As a function of r and s, R is a solution of the differential equation adjoint to (A): 


R,, — | (2) + ( = ) |- 0. 
r+s/, Jy a « 


(b) Onr = é, R, — AR/(E + 8) = O, and 
ons = 7, R, — A\R/(r + n) = 0. 

(c) R(é, n; &, n) = 1. 

It is true in general that if R(r, s; £, 7) is the Riemann function for a given differential 
equation, then R(é, 7: r, s) is the Riemann function for its adjoint. The Riemann func- 
tion for the differential equation adjoint to (A) was obtained by Riemann and is given in 
[1, Sec. 82]. Interchanging parameter and variable points in that solution gives the 
Riemann function wanted here. The result may be checked of course by substitution 
into (a), (b), and (c) above. 


; | 
ir. a8. «) = rte) | - ; , Be 9 
Hr, 83, ez MIM MG -Et e+e d’ 


where F is the hypergeometric function. If \ is integral, F reduces to a polynomial. 
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It is shown in [1, Sec. 72], that the entropy change through a shock wave is of the 
third order in the shock strength. Thus shocks here must be weak enough that third 
and higher order terms are negligible. The velocity of a shock wave is then 


y¥+1 (y + 1)? (u — w)’ 
4 32 Co , 








U =u +e + (u — U) + 
where the upper signs are valid for a forward facing shock in which U > up , the lower 
signs for L’ < up» , zero subscripts refer to conditions ahead of the shock, and no sub- 
script to the conditions behind. 

III. The transmitted wave. The two fluids are at rest initially except that a wave 
of known properties in the left-hand fluid, say, is moving toward the interface. 

The initial state of the right-hand fluid may be represented by a point, A, in the 
(l’, wu) hodograph plane. (Primes are used for quantities in the right-hand fluid, no mark 
for those in the left; e.g. l’ is the 2 function for the right-hand fluid.) If the transition 
from rest to motion is continuous, characteristics crossing the transition, the T— 
characteristics here, must all go through A. Thus the whole flow of the right-hand fluid 
maps itself onto the !'— characteristic through A and so is a simple wave. The assump- 
tion of weak shocks implies that the simple wave flow is also correct in the discontinuous 
shock case [1, Sec. 72]. Letting sf be the constant value of s’ in the simple wave, one has 

l’ —u = 2si (1) 
throughout the right-hand fluid. 

The solution process carried out below determines (¢, u, and l’ consistent with Eq. 
(1) along the interface 7. More definitely, ¢ and u are determined as functions of a con- 
venient parameter q along J. Then, since / is a particle path, x = x(q) is found by inte- 
grating 

dx dt 


= uz 


dq dq 
The C+ characteristics in the simple wave satisfy 
dx = (ett r+r53 x) dt, 


where different constant values of r’ give the different characteristics. This equation 
may be integrated under the condition that z = x(q), ¢ = ¢(q), andr’ = r’(q) = (1/2) 
[l’(q) + u(q)] where the characteristic crosses J]. The integration gives 


x — x(q) = Ez r’(q) + z = Ba te — t(q)], (2) 








where zx and ¢ designate an arbitrary point on the C+ characteristic. 

If the flow is continuous it is completely described by Eq. (2). At any point (2, ¢), 
Eq. (2) determines the corresponding value of g. This gives r’ = r’(q). Also s’ is known 
to be the constant sf. Then u = r’ — s/ andl’ = r’ + s{ so that the flow is determined 
at that point. 

The flow resulting from a sharp drop in pressure at initial point B, where the front 
of the oncoming wave in the left-hand fluid first reaches the interface, is also described 
by Eq. (2). The pressure drop means a corresponding range of values of r’ is concentrated 
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at B. This fits the algebra of (2) when a range of values of parameter g is assigned for 
which ¢(q) = ¢, , the initial time, but r’(q) is variable and covers the necessary range in 7’. 
The flow resulting from a sharp rise in pressure at B is discontinuous with a trans- 

mitted shock front preceding the simple wave region described by Eq. (2). In this case 
it remains to determine S, , the path of the shock front transmitted into the right-hand 
fluid. This path, x = X(q) and ¢ = T(q), is found by integrating 

dX . aT ; 

dq = T dq ’ (3) 
where U; is the transmitted shock front speed. U7 is a known function of flow speeds 
and pressures on the two sides of the shock and so is a known function of g here. X and 
T are related by Eq. (2) with x = X andt = T. Eliminating X between Eqs. (2) and (3) 
then gives a first order ordinary differential equation for 7 = T(q). 


, aT adx(q) yy’ +1, vy’ —3 ,|| dT dt(q) vy’ +1 aw dr’(q) 
—_—— — = | —— 7'( —_—— ¢§ — — ——— | + —_—_—_ — t(q)| —~- 
Ur dq dq | i Q + 2 ' |Ldq dq bi 2 U7 he dq 





Solving this gives 7 (¢). Equation (2) then gives X(q) and so the shock front path S; . 

The above solution of the transmitted wave problem is valid so long as character- 
istics of the same kind do not intersect. For at such an intersection two distinct values 
of r’ are given and so the flow is ambiguous. If the oncoming wave is such that the 
pressure along the interface decreases as time goes on, the intersection of C+ character- 
istics does not occur. For Eq. (2) may be written 

, 
xz — x(q) = a [l’(q) — 2Qs{][t — t(q)]. 
Since 1’ decreases as the pressure decreases, the slope of these straight lines increases 
with time and so they cannot intersect to the right of the interface in the (7, ¢) plane. 

On the other hand, if the pressure increases along any section of the interface, the 
C+ characteristics through that section do eventually intersect. The flow ambiguity 
is eliminated by the formation of shock discontinuities. It is assumed here that the 
region of interest does not contain such intersections. 

Since the flow in the simple wave is determined by the C+ relation, a single value 
for the slope of the C— characteristic is obtained at each point. Only one C — character- 
istic may then pass through each point. The C— characteristics do not intersect. 

IV. The interface. The interface lies along the I characteristic through 1 in the 
(l’, x) hodograph plane. Now I’ is a monotonically increasing function of p in the right- 
hand fluid with p in turn a similar function of 1 in the left-hand fluid. Thus, along J, 
l’ is a function of J, l’ = f(l). Equation (1) valid throughout the right-hand fluid then 
becomes 

f(l) —u = Qs{, 


the equation of the J curve in the (I, u) hodograph plane of the left-hand fluid. This one 
interface relation is equivalent to a pair of relations 
u=u(q) and l= (U(q) 


or 
r=r(q) and s = s(q) 


for any convenient choice of parameter 4. 
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The boundary condition along J is that the interface is a particle path: 


dx _ dt 


dq = Uu dq 


This may be expanded in characteristic coordinates to 


dr ds dr ds 
p — aes P — zs rc — § — -t- — . 
-" dq — ve dq v al « dq t, ds] 


Eliminating x by relations C+ and C— finally gives 


14242 we, 
dq dq 

If B is the point at which the oncoming wave first meets the interface, the C+ 
characteristic through B is the front edge of the wave. This characteristic continues 
ahead through the vertex of a reflected rarefaction fan or through a shock front, as 
the case may be, on to the interface. At B; , i.e., at B on the interface, conditions are 
those along the C+ characteristic through B. In the hodograph plane this means that 
B, lies at the intersection of the [+ characteristic through B and the J curve. 

The relative positions of A, B, and B, determine the local nature of the discontinuities 
in both fluids at B, for the jumps from A and B to B, in the hodograph plane represent 
those discontinuities. If B, is to the right of A (as sketched in Fig. 1), 1 and therefore 
the pressure rise in going from A to B,; ; a shock front has been transmitted into the 
right-hand fluid. If B, is to the left of A, a centered rarefaction wave has been trans- 
mitted. Similarly, if B; is to the right of B (as in Fig. 1), a shock wave is reflected into 
the left-hand fluid. If B,; is to the left of B a rarefaction is reflected. The two types of 
reflection are discussed separately in the two following sections. 











Fic. 1. Hodograph for left hand fluid. 


For the case of a transmitted shock wave, with B, to the right of A as in Fig. 1, 
D. C. Pack [6] has derived convenient analytic expressions for determining whether a 
shock or a rarefaction is reflected back into the left-hand fluid. One may rederive these 
relations geometrically as follows. 

For a reflected shock, from Fig. 1, 1(B) < 1(B;). Since the pressure is an increasing 
function of 1, p(B) < p (B;) and p(B) — p(A) < p(B;) — p(A). Since, also from Figure 
1,0 < u(B,) < u(B), the pressure difference inequality is strengthened in the form 

p(B) — p(A) _ p(Br — plA). 


— u(B) u(B;) 
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The Rankine-Hugoniot relations for a shock of velocity U, where the subscripts 1 
and 2 refer to the two sides of the shock, read 


p(U — u,) = p{U — wu) 
D2 — Di = p(U — uw)? — p(U — uw)? = p(U — w)(u2 — uw). 
Specializing the last relation to the present situation gives 
p(B) — p(A) = poUu(B), 
p(By) — p(A) = piUu(B), 
where pp, and pj, are the densities of the undisturbed fluids and U, and Uy are velocities 


of the oncoming and transmitted shocks, respectively. Substitution gives Pack’s first 
inequality: For a reflected shock 
poo < por . 
Pack’s second inequality refers to a strong oncoming detonation wave front for 
which p(A) is negligible with respect to p(B) and the detonation velocity, D, satisfies 


the Chapman-Jouguet condition D = u(B) + c(B). These assumptions together with 
the Rankine-Hugoniot equations above and the speed of sound relation, c’ = yp/p, 


imply 
D = (y + 1)u(B), 
p(B) = po Du(B), 


II 


and 
iB) (4 - 1) =u, 
Po py 


equations needed below. Here u, and p} are the particle velocity and density behind a 
fictitious shock front of pressure p(B) propagated into the undisturbed right-hand fluid. 

For the reflected wave to be a shock, Fig. 1 shows that point B lies above the shock 
polar of the right-hand fluid, the locus of points to which a shock could move point A. 
In particular B is above the point on the shock polar with 1 = /(B). In symbols 


u(B) > u;. 
Since u(B) and u, are both positive one has 
[u(B))? > uz 
and, with the equations listed above, this last may be reduced to Pack’s second in- 
equality: For a reflected shock, 
l ] 
(fy + Dalla — 3) < 1. 
Po ps 
V. The reflected rarefaction case. If B is to the right of B, in the (1, u) hodograph 
plane, a centered rarefaction wave is reflected back into the left-hand fluid at the inter- 


face (see Fig. 2). The front of the oncoming wave, the C+ characteristic through B, 
goes on through the vertex of the centered rarefaction to B,; . The value of ¢ is then 
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known along BB, in either the hodograph or the characteristic plane, this value being 
just the constant { = t at B. If BD is the last C— characteristic through B unaffected 
by the reflection, the value of ¢ along BD is also known. (In the hodograph and character- 
istic planes D lies somewhere on the C— characteristic through B. It is shown to the 
left of B in the hodograph sketch, Fig. 2(b), for definiteness.) 


(a) | c- 
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Fic. 2. Reflected rarefaction: 
(a) Physical plane 
(b) Hodograph 
(c) Characteristic plane 


The value of ¢, i.e., the solution of partial differential equation (A), at any point P 
in the characteristic rectangle determined by B, B, , and D is then given following 
Reimann [2, Chap. 5, Sec. 4] by 
“(P) = (B)R(B; P) 4 [ RI; P\(t + — ds + “RU; P(t + = dr, (4) 

‘ dati = rae ee © "ee " * 
where / and F are the points of intersection of the characteristics through P and BB, 
and BD, respectively, J is a general notation for the moving point of integration, and 
R is the Riemann function. In particular ¢ may be determined along B,G, the last 
characteristic in the centered rarefaction wave. 

Another characteristic rectangle, B;GJK, is determined by B,G and interface curve 
I considered as a curved diagonal. (Since a characteristic such as B;G may act as a 
“transition” line, a line along which the map of the physical plane into the hodograph 
plane is folded, there is no difficulty if J cuts back into rectangle B,;,GDB. It is only 
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necessary to treat separately regions where the slope of J is greater than or less than the 
slope of the !'— characteristic.) 

Values of ¢ at any point P in rectangle B,GJK are determined as before in terms of 
values of ¢ along B,G, where they are now known, and along B;K, a characteristic in a 
fictitious region of the left-hand fluid continued across the interface. If E’ and F’ are the 
intersections of characteristics through P with B;K and B,G, respectively, then 


FP’ 
RJ; P)(, + ; Mt) dr. 


toe wane ai At 
t((P) = t(B)R(B, ;P) + RJ; p(s, + ) ds + | a 


/ Br / By 


The partial derivative of ¢ with respect to r is obtained by computing the change in ¢ 
from P to nearby point P’ on the same s = constant line. 





t,(P) = 1(B)RAB, ;P) + RF’; Py(1, + Mt \ ey 


ie ae 
E’ } t nF’ ; ; ; 
+ | RJ; P(t pay | R(J;P)(t +- a di 
“BI ie J Br is aa 
where parameter point P has coordinates (é, ») and 
> OR 
“a =~. | 
: Ge le «% 
Similarly 
7 w 1. 
t,(P) = (B)R,(B; ;P) + RE P(t ae Mt \ Kg ) 
r+s 
E’ Fe 
ie At Ne 
i RJ; eer I RAJ; P)\t rarer Ro 


Considering P to lie on the interface and substituting these forms into the interface 
boundary condition, 


, dr ls 
I: t.(P) = (P) — t,(P) = (P) = 0, 
dq dq 
gives an equation to be solved for the unknown values of ¢ along B,K, that is the value 
of ¢ at variable point E’. The equation is a first order integro-differential equation in 
one independent variable. For, the interface lies on the known curve 


r=r(q), 8 = s(q). 


Point P has this r(q) and s(q) for coordinates. Fixed point B,; has known coordinates 
rr and s; . Then F’ has coordinates r(q) and s, , while E’ has coordinates r; and s(q). 
The coordinates of all points and so the value of all functions thus depend on known 
constants and the one variable g. A step-by-step finite difference numerical process in 
one dimension can then be used to approximate the solution of the equation. 

With ¢ so determined along B,;K, Riemann’s formula gives ¢ at any point P within 
characteristic rectangle GB,;K J. Thus ¢ is known as a function of r and s in the reflected 
wave. From the values of ¢ along J, ¢ is known throughout the transmitted wave as was 
pointed out previously. Corresponding values of x are found by integrating the character- 
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istic relations in either fluid. To obtain the path of the interface, the boundary condition 


is integrated to give x = 2(q) knowing ¢ = ¢(q). 

VI. The reflected shock case. If B is to the left of B; in the (l, u) hodograph plane, 
a shock discontinuity is reflected back into the left-hand fluid at the interface. The 
shock curve Sz in the (2, ¢) physical plane becomes two curves S, and S, in the (J, ~) 
hodograph plane (see Fig. 3). S, gives conditions just ahead of the shock and so starts 
at B. S, gives conditions just behind the shock and so starts at B, . 
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Fic. 3. Reflected shock: 
(a) Physical plane 

(b) Hodograph 

(c) Characteristic plane 


The homentropic assumption implies a weak shock across which the C+ character- 
istics are continuous. Thus r is the same on both sides of the shock and is a convenient 
parameter for the shock curve. The shock boundary condition then reads 


where U, , the velocity of the reflected shock, is a function of flow speeds and pressures 
on both sides of the shock. If s, is the value of s in the oncoming wave just ahead of the 








24 C. M. ABLOW (Vol. XVIII, No. 1 


shock, for a given r, and s, the corresponding value of s just behind the shock, Uz is 
the function of r, s, , and s, given in Sec. II above. 

The oncoming wave flow is “‘known’’, i.e., x and ¢ are known functions of r and s. 
The shock boundary condition becomes a differential equation for s, : 


S, $2, + 2, ds, = Ud(t + ¢t, ae.) 3 
dr dr 


where the partial derivatives z, , etc., are all evaluated at (r, s,). 
Applying the shock condition just behind the front gives 


) ds, _ Al des), 
%, + 2, _* Ur\t, + t, a 


The partial derivatives are related by the C+ and C— characteristic conditions so that 
x, and x, may be eliminated to give 
' <—s y¥+1 : + y¥—3 ds, 

S, (1 2+ r+ vt a)e + (« e~ eee = Sra 2 e. 
Here #, and ¢, are to be evaluated at (r, s,). However, the forms of ¢, and ¢, as functions 
of r and s, are not yet known. 

Now ¢ may be viewed as a known function of r and s, or as an unknown function of 
rand s, . Equating these values of ¢ along the shock then gives 


S, : r,s.) = Ur, %). 


The problem restated is then to solve the partial differential equation (A) for ¢ asa 
function of r and s in the angular region between 7 and S, in the (r, s) characteristic 
plane subject to conditions J along the interface and S, , S, , S; along the shock (see 
Fig. 3(c)). 

As was done with the reflected rarefaction, it is convenient to consider the flow 
continued outside the angular region between J and S, to the characteristics through 
B,. Forms for ¢, ¢, , and ¢, given in the previous section may be substituted in the various 
boundary conditions as before. 

The substitution into interface condition J gives an integro-differential equation 
connecting ¢ at F and ¢ at E, . Substitution into S, connects ¢ at F with ¢ at £, . But in 
the step-by-step process by which the equations may be solved, ¢ will already be known 
at EL, so that the S, equation actually relates ¢ at F and function s, = s,(r). (The r co- 
ordinate of F, r, is now the parameter q in terms of which all points are expressed.) 
Finally, substitution into S; gives an equation relating ¢ at F, s,(r), and s,(r). The four 
equations J, S, , S, , and S; thus give the necessary number of relations to determine ¢ 
along B,F, ¢t along B;E, s,(r), and s,(r). 

As in the case of the reflected rarefaction, knowledge of ¢ along characteristics B,F 
and B,E and of the location of S, and S, in the (r, s) plane permits the determination of 
t everywhere and so, using either the characteristic relations or the appropriate boundary 
condition, the corresponding values of x everywhere. 

The next section contains an application of the present method to a specific example. 

VII. Refraction of a detonation. A problem of some interest is the determination 
of the flow of a material due to the detonation of a layer of high explosive on its surface. 
An idealized form of such a problem, amenable to the present methods, reads as follows. 

A semi-infinite slab of material filling the space x > 1 is given. A layer of high ex- 
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plosive of unit thickness covers the slab. At time ¢ = 0 the explosive is detonated every- 
where on its free surface (x = 0) simultaneously. The gaseous combustion products are 
assumed to expand into a vacuum. 

The detonation front burns into the high explosive at known speed D. It is assumed 
that the detonation is a Chapman-Jouguet process, i.e., that particle and sound speeds 
just behind the detonation are related by 

ute= D. 


This and the equations for the conservation of mass and momentum determine all the 
flow quantities immediately behind the detonation front in terms of D and the known 
initial density of the high explosive. Assuming the combustion products form a polytropic 
gas with zero density at zero pressure fixes a, 8, and y, the constants in its equation of 
state. 

The slab material is also assumed to be a polytropic gas with known equation of 
state. Some discussion of the experimental basis for this last assumption is given in 
[4, 6, and 7]. 

Since along the detonation front 


dx = Ddt 


(u + c) dt, 


this front lies on a C+ characteristic. The flow field behind the front is covered by 
C— characteristics. The detonation front carries constant values of pressure and flow 
velocity and so is represented by the single point B in the (1, uw) hodograph plane. The 
field behind it must then all lie on the [— characteristic through B and so is a simple 
wave. In the wave s = 8 , a constant. 

In the physical (z, ¢) plane the whole wave starts at the origin, z = ¢ = 0, and soisa 
centered simple wave. Straight lines through the origin are all C+ characteristics and 
correspond to points on the ['— characteristic between B and V, the (vacuum) point 


where p = 1 = 0. 
The C+ characteristic equations may be integrated here to give 


r= (141,475 3,), (5) 


— 





This together with 

8 = &% 
describes the centered simple wave in the burned gas, the oncoming wave of the preceding 
theory. 

Since a simple wave is a degenerate form, some simplifications of the reflected wave 
boundary conditions occur here. In this example it is assumed that B lies above inter- 
face curve I in the (l, u) hodograph plane so that the reflected shock considerations of 
Sec. VI apply. Points just ahead of this shock lie in the simple wave. The S, curve in 
the hodograph plane thus coincides with ['— characteristic through B. Boundary con- 
dition S, is no longer needed to locate S, . When the other conditions have been used 
to determine ¢ as a function of r along the front, S, , or its equivalent 


dx , dt 
ar _— Ura» (6) 
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may be used to find x as a function of r and so the location of the reflected shock in the 


(x, ¢) physical plane. 
Boundary condition S, is needed as before to locate S, , the hodograph curve of 


points just behind the reflected shock. 
One may eliminate x between the shock condition, Eq. (6), and the equation valid 


in the oncoming simple wave, Eq. (5), to obtain 


r+1 (2+4 1~3 ») _ yy at. 
. + - 8 “oe ee Mee 


By boundary condition S; , the values of ¢ at corresponding points of S, and S, are the 
same. Hence ¢ and its derivatives in the equation may be considered to lie along S, . 
The shock boundary condition then reads 
; y¥+1 . y+ 1 +<3,) ds, 
S Te - == ——_—-— — PP == -« ——s © = — 
i. 5 t Us; 9 T 5 t T S de ’ 


where /, ¢, , and ¢, are evaluated along S, . 
Finally the interface boundary condition, with parameter q identified with character- 
istic coordinate r of the left-hand fluid, reads 


ds, 


ee ae - 
r 8 dr 


The forms for the partial derivatives of ¢ to be substituted in S, , S, , and J may be 


conveniently written 


Il 


t,(P) = t,(F)R(F; P) + ®,(E;P), 
t,(P) = ti(E)R(E; P) + #,(E; P), 





where 
BF) = ¢(F) = = x 
. 4 = 7 = - r. Sr) 
F\ ) r / dr r, 87), 
iE) = t(£) = £ Ce) 
6,(E: P) = AF) R(F; P) + «(B)RAB, ; P) 
r+ 8; 
. At i: \t 
RAJ; P (1, —— | ds RAJ; P)(1, + —“ -) lr, 
+f . ’ TT +a + | = ) rt sg, o 
and 
. AL(E) : meee 
@,(E£;P) = —— RE;P) + (BRB; ; P) 
rrt+s 


! [Rus P ( ee le + [RQ Pyle jm ) | 

+ dif ht, + ————] @& J; > ———) &. 
JBr > rr+s JBr =. r+ 8, 

In using these formulas, P is replaced by either Ps or P; , points with the same r co- 

ordinate as F and lying on S, or J (see Fig. 3(c)). Similarly # is replaced by Fs; or E, , 

points whose r coordinate is r; and whose s coordinate agrees with Ps; or P;. 
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Making the substitutions into S, , S,, and J and then solving for the derivatives of 


the unknown functions gives 


U(F) = [®, — (1 + &)®(Es ;P3)/( + &)R(F; Ps) (7) 
—— inieaieaiiaiatie te . . ds;| . ds; : 
tp(E)) = tp(F \RF } P;) + ®,(E, + 1 _ $,(E, ; Py) yy ~~ - RE; ; P;) (8) 
& (P,) = 4,4,/(1 + &)[(EIREs ;Ps) + (Es ;Ps)], (9) 
where 
y+1 , +1 r= 3,,) 
se = eS eae ee ak 
®, . uP) / (Us 9 r 9 So 
7- + 1 ? 1 — 3 
” (ue + 1537474 44)/(u-4 r— 75 s) 
and 


: y¥—1 2 &—-& 
Ui (ioe Soe — +49) 1 2 (2 =) 
_&t+ v" ( 2 & =*)"] 
32. Wy —iIrts/ J 
Equations (7), (8), and (9) then constitute three simultaneous, ordinary, first order, 
integro-differential equations for tp , tg , and s, , the values of ¢ along characteristics 
B,F and B,E and the location in the hodograph plane of points just behind the reflected 
shock. Using these, the value of ¢ at any point of the reflected wave may be computed. 
It is necessary to calculate ¢ along s, using Eq. (4) as the integration proceeds for sub- 
stitution into ®, . The computation of ¢ along J may be carried out when convenient. 
Appropriate integrations then give x along s, or J and so the path of the reflected shock 
or interface in the (x, t) plane. The path of the transmitted shock may then be computed 
as described in a previous section. 
VIII. A numerical example. The theory of Sec. VII has been applied to determine 


the constants in an equation of state for the gaseous combustion products of a high 
explosive known as Composition B (See [7].) The equation is assumed to have the form 





p = ap’ — B. 
At the detonation front one has the Chapman-Jouguet condition, 
uw t+aQ = D, 
the conditions for conservation of mass and momentum, 
plu. — D) = —p D, 


and 
Pi = po Du, , 


and the definition of the speed of sound, 


¢, = (ayp)')"”. 








28 C. M. ABLOW [Vol. XVIII, No. 1 


These relations may be combined to give 


B = [p) D? — (y + Ipil/y 


and 
a = (po D*/y)[y(8 + po D*)/(y + po D*)""". 


Thus two of the constants in the equation of state, a and 8, are determined in terms of 
the third, y, the initial density of the high explosive, p) , and the known pressure and 
velocity of the detonation, p, and D. 

For various values of y, curves of U; , the transmitted shock velocity, versus (x — 1), 
the distance into the slab, may be plotted as in Fig. 4. These curves may then be compared 
with the experimental points given in Table I of Ref. [7]. The curve for y = 11/9 is 
seen to agree fairly well with the experimental points near the surface of the slab. This 
value of y is also suggested for combustion products in Ref. [1], Sec. 38. A more complete 
comparison of theory with experiment has been made by W. E. Drummond [10]. 
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Fia. 4. Transmitted shock velocity vs. depth. 
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As y varies, the (p, p) curves all pass through the same point (p; , p,) with the same 
slope, c; . The equation of state has been taken in polytropic form to permit use of the 
Riemannian theory. For 8 ¥ 0 the equation can only be an approximation valid in the 
neighborhood of p = p; , p = p; . Some indication of how large this neighborhood may 
be is given by (0), the density for which p = 0, since this density would be zero in a 
correct equation of state. 

The equation of state of the metal slab is also in polytropic form. This form has been 
given some theoretical standing by Murnaghan in [5]. Constants used in the equation 
are those for aluminum given in [4]. 

The table below lists numerical values of interest. Densities are given in gm/cm’, 
pressures in megabars, and velocities in cm/microsec. 











po = 1.67 pr = 0.26 y’ = 4.267 
D = 0.77 pi = 2.263 6’ = 0.171 
c = 0.567 a’ = 0.00214 
Y B a p(0) 
2.81 0 0.0263 0 
2 0.104 0.0715 1.28 
5/3 0.177 0.1124 1.31 
7/5 0.261 0.168 1.37 
9/7 0.307 0.201 1.39 
11/9 0.336 0.222 1.41 





The calculations have been ably carried through by T. C. Poulter, Jr., on an internally 
programmed digital computer. The work has been supported in part by the U. S. Army 
Bureau of Ordnance. 
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BOOK REVIEWS 


(Continued from p. 14) 


grated form. The method is then applied to the determination of the eigen-values for the differential 
equation of the vibrating string. 

Among the most interesting and useful parts of the book are the collections of exercises and research 
problems at the end of each chapter which give an indication of the truly enormous range of potential 
applications, including such disparate items as a comparison of the arithmetic and geometric mean, 
Fibonaccian search, the design of missiles and the crossing of a desert by jeep. 

The presentation is keyed to the level of ordinary calculus and is kept informal to the point where 
important proofs are left to the reader or are deferred to a second volume which is to be on a more 
technical level. These further developments will be looked forward to with keen expectation. Mean- 
while it is a pleasure to recommend to anyone whose job is the application of mathematical tools in 
any context whatsoever the study of this most stimulating and rewarding book. 


Martin J. BECKMANN 


Theory of relativity. By W. Pauli. Translated from the German by G. Field. Pergamon 
Press, New York, London, Paris, Los Angeles, 1958. xiv + 241 pp. $6.00. 


Pauli’s ‘‘Theory of Relativity”’ is a remarkable book. It was written in 1920 when the general 
theory of relativity, the special theory of relativity and the author were respectively four, fifteen, and 
twenty-one years old. Now, almost thirty years later, it is still one of the best books on the subject. 
The point of view remains fresh; Pauli’s comments and criticisms almost invariably remain relevant. 

One of the striking features of the book is the description of the historical developments preceding 
the formulations of the special and general theories of relativity, the discussion of rival theories and of 
the relevant experiments. Pauli’s book covers the field of relativity theory more extensively than most 
texts. Although written in an encyclopedic style which omits many derivations and proofs, enough 
details are given to enable the mature reader to supply the missing steps, with only occasional use of 
some of the many literature references. 

“Relativitaitstheorie’’ was written as an article for the Encyklopidie der Mathematischen Wissen- 
schaften (Vol. 5, Part 2; No. 19), appearing there in 1921, and shortly afterwards was reprinted as a 
book. The present volume consists of an excellent translation by G. Field of the original text of “‘Rela- 
tivitdtstheorie,’”’ totalling 206 pages, and of 26 pages of supplementary notes in which Pauli reviews 
selected developments in relativity theory which occurred from 1920 to about 1955. 

The general theory of relativity has had a peculiar position in modern theoretical physics. Gravita- 
tional forces are extremely weak. They can be ignored in most of atomic and nuclear physics and, on 
the astronomical scale, Newtonian gravitational theory, which is not even Lorentz invariant, agrees 
so well with observation that the few minute deviations predicted by general relativity are very difficult 
to measure. As a result, physicists have paid little attention to general relativity between 1920 and 
1955, and this is reflected by the fact that Pauli can do reasonable justice to developments during this 
period in some 20 pages. 

Since the mid-fifties, there has been renewed interest in gravitation. This is partly due to a general 
interest in field theories, spurred on by the successes and limitations of quantumelectrodynamics, and 
partly to the possibilities of achieving greater aceuracy in astronomical observation when telescopes 
can be lifted above our atmosphere or when suitably chosen motions of artificial celestial bodies can 
be studied. The resurgence of interest in general relativity theory is reflected by the fact that in the 
few cases where Pauli’s information is outdated, this is due to very recent work, such as that of Bass 
and Pirani (Phil. Mag. 46, 850; 1955) on the relativity of centrifugal force, and that of Robinson and 
Bondi on gravitational waves (Nature 179, 1072; 1957). 

ALFRED SCHILD 


(Continued on p. 36) 
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THE PHYSICAL REALIZABILITY AND REALIZATION OF LINEAR 
PHASE SHIFT NETWORKS* 


BY 
PAUL M. CHIRLIAN 
New York University 


Abstract. Necessary and sufficient conditions for the physical realizability of linear 
phase shift networks are imposed upon the network transfer function. A technique is 
presented for approximating non-physically realizable transfer functions with those 
which are physically realizable. It is demonstrated that if the approximation is properly 
limited, the unit impulse response can be made non-negative for all values of time. 
These results are extended to networks whose transfer functions phase shift is the sum 
of two terms, one which varies linearly with frequency and a second which varies in 
discrete steps of m radians. 

Introduction. Investigations of the transient response of linear networks have led 
to the conclusion that a network transfer function whose phase shift is a linear function 
of frequency is often desirable. Of course, if such a network is to be built, it must be 
physically realizable. Necessary and sufficient conditions for this physical realizability 
have been obtained. The effect of these conditions on the transfer function can be readily 
visualized. 

If a desired transfer function is not physically realizable, then it should be approxi- 
mated by a transfer function which is physically realizable. An approximation of this 
type is presented, and it is shown that the accuracy of the approximation is improved as 
the magnitude of the slope of the phase shift (the network time delay) is increased. 
The unit impulse response of the approximating network can always be made a non- 
negative function of time if the given transfer function is of finite bandwidth and the 
approximation is properly limited. 

The following definition will be used throughout the body of the paper. The ratio, 
expressed as a complex function of frequency, of the steady state response of a network 
to its sinusoidal input will be called the transfer function. It will be written as T'(w) 
exp [j@(w)], where 7T'(w) is the amplitude function and @(w) is the phase function. If 
#(w) = — kw, where k is a real positive constant and 7'(w) is a real positive function of 
w, then this is the transfer function of a linear phase shift network. If 7J(w) takes on 
negative as well as positive values, then this will be called the transfer function of a 
stepped linear phase shift network. The term stepped is used since the negative values 
of T(w) may be accounted for by adding discrete steps of z radians to the phase function. 
The expression, generalized linear phase shift network, will be applied to networks 
which are either linear phase shift or stepped linear phase shift. 

It shall be assumed that any amplitude functions considered here satisfy the Dirichlet 
conditions and, hence, possess a Fourier transform. 

The physical realizability of generalized linear phase shift networks. If an arbitrary 


*Received December 8, 1958. This paper is based upon a portion of a thesis which has been accepted 
by the faculty of the Graduate Division, College of Engineering, New York University, in partial ful- 
fillment of the requirements for the degree of Doctor of Engineering Science. 
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transfer function is assigned to a network, it is often found that the response of the 
network to a suddenly applied signal must appear before the signal is applied. Such a 
network is not physically realizable. In fact, the definition for physicai realizability 
which shall be used in this paper is: A network is physically realizable if its response 
to a signal, applied at time ¢ = 0, is zero for all ¢ < 0. Hence, using the convolution 
integral, a necessary and sufficient condition for physical realizability is that the response 
of a network to a unit impulse applied at time ¢ = 0 is identically zero for all values of 
time less than zero. In order to convert this condition into one which can be interpreted 
physically, some properties of generalized linear phase shift networks will be discussed. 
The unit impulse response W(t) of a generalized linear phase shift network can be 


written as 


, tw. 
Wi) =- | Tlw) cosw(t — &) de. (1) 
T Jo 


It can be seen that W(t + k) is an even function of time. Hence, if a generalized linear 
phase shift network is to be physically realizable, it is necessary and sufficient that 


W(it+k) =0 for t< —k and i> k. (2) 


It is of interest to note that the unit step response of a physically realizable generalized 
linear phase shift network is such that 
, _ [7T0), t>k 
A(t +k) =} 


0, t< —k 


From the Fourier inversion relation we obtain 


Ita) = 2 | W(u + k) cos wu du, (3) 
where u = ¢ — k. Thus, we may state the following theorem: 
Theorem 1. A necessary and sufficient condition that a generalized linear phase shift 
network, whose transfer function is T(w) exp (— jkw), be physically realizable is that T(w) 
be expressible by the following Fourier transform: 


. 


Tw) = | H(u) cos wu du. (4) 


“0 

Proof. The necessary condition is proven by assuming that W(t) = 0 for ¢ < 0 and 
then applying Eq. (2) to Eq. (3) and then substituting H(u) = 2W(u + k). The suf- 
ficiency condition is proven by hypothesizing Eq. (4), substituting H(u) = 2W(u + k) 
and then comparing it with Eq. (3). Hence, W(u + k) = Oforu > k, but W(u + k) 
is an even function of u. Therefore, W(t) = O fort < 0. 

For an interpretation of Theorem 1 we can make use of the sampling theorem of 
Shannon [1]. From this theorem it can be seen that 7'(w) is completely specified by a 
set of points spaced 2/k radians apart. Hence, increasing the magnitude of the phase 
slope, k, increases the accuracy of the approximation, in the sense that the function 
may be specified at more closely spaced points. It should be noted that the constant 
k is also the delay time of the network. Thus, the maximum value of k may be limited 
by the design constraints of the system. 
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An approximation procedure. When an arbitrary generalized linear phase network 
transfer function 7',(w) exp (— jkw) is specified, often 7,(w) does not satisfy the con- 
ditions imposed by Theorem 1. A technique is given here for approximating the given 
transfer function with a physically realizable transfer function, T.(w) exp (— jkw). 
The form of the physically realizable approximate amplitude function, T,(w), will 
depend upon the criterion of approximation used. In this case the Fourier criterion was 
used. That is, if 7',(w) is of the form, 


T(w) = i H,(u) cos wu du, (5) 
0 
where 
— 
Hw) = = [ T(s) cos wu deo, (6) 
then 7',(w) will be defined as 
T.(w) = [ H,(u) cos ts du. (7) 


An increase in k will, of course, increase the accuracy of the approximation. Note 
that an increase in k will not just introduce a constant time delay in this case because 
T ..(w) will also vary. 

Approximation procedures which guarantee a non-negative unit impluse response. 
It is often desirable for the unit impulse response of a network to be non-negative for all 
values of time (a monotonically increasing unit step response). If a transfer function 
of a linear phase shift network, 7,(w) exp (— jkw), which does not have a non-negative 
unit inpulse response, is specified, then often it is desirable to approximate T,(w) with 
a physically realizable approximate T’,(w), such that the resulting network has the desired 
non-negative unit impulse response. Procedures for obtaining such a 7,(w) will be 
presented. Two of these require that 7',(w) be bandwidth limited. That is T,(w) = 0 if 
w > w, , where w, is some specified angular frequency. In many instances such low pass 
characteristics are specified or closely approximate the specified characteristics. 

In the following theorems it is assumed that the networks are generalized linear 
phase shift with a phase slope of magnitude k. 

Theorem 2. If T,(w) is a physically realizable approximate of some T(w) where T,(w) = 
Oifw > w,,7,(w) > 0 for0 < w < w, , andw.k < 2/2 then the unit impulse response 
will be non-negative for all values of time. 

Proof. From Eq. (6) we obtain 


ie : [ " Filed eneensidhe. 
0 


For the range of integration, 7,(w) > 0 and cos wu > 0, since 0 < wu < wk < 4/2. 
Hence, H,(u) > 0 for all | u | < k; but (from the proof of Theorem 1) W(t) = H(t — k)/2 
for 0 < i < 2k and W(t) = Ofort < Oandt > 2k. Therefore W(t) > 0 for all ¢. 

Thus, if the magnitude of the phase slope is equal to or less than 2/(2w,), the unit 
impulse response will be non-negative. It may be desired to obtain a larger phase slope 
magnitude, without altering the shape of the transient response (although time delay 
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will be introduced). This may be accomplished by not increasing the upper limit of inte- 
gration of Eq. (7) so that now T’,,(w) is given by 


w/(2wc) 


Tw) = [ H,(u) cos wu du. 
v0 

If, in addition to the requirements imposed in Theorem 2, the specified amplitude 
function 7',(w) is a monotonically decreasing function of time, then the restrictions on 
k imposed in Theorem 2, can be relaxed somewhat. 

Theorem 3. If T.(w) is a physically realizable approximate of some T,(w) where T,(w) = 
Oifw > w,,7Ti(w) = Ofor0 < w < w, , T:(w) ts a monotonically decreasing function of 
w, and w.k < x then the unit impulse response will be non-negative for all values of time. 

Proof. The proof follows that of Theorem 2 except the fact that 7',(w) is a mono- 
tonically decreasing function of w allows cos wu to take some negative values, provided 
wu <7. 

Thus, a better approximation to 7, (w) can be obtained when 7’, (w) is a monotonically 
decreasing function of w. 

If the specified amplitude function 7,(w) is such that 


© 


| T,@) coswidw > 0 forall O<u<M, (8) 
#0 

then any physically realizable approximate of 7',(w) will result in a non-negative unit 

impulse response provided that k < M. An example of this is a 7',(w) which decreases 

steadily to zero as w approaches infinity and is convex downwards. The integral of 

Kq. (8) will be positive for all values of u, (see Titchmarsh [2]) Thus, any positive value 

of k may be used. 

As far as the problem of actually realizing these networks is concerned, it can be 
shown that these networks cannot be realized with a finite number of lumped linear 
elements. However, Corrington and Sonnenfeldt [3] and Kallman [4] have developed 
methods for realizing linear phase shift networks, using distributed elements. These 
techniques can be applied to generalized linear phase shift networks. 

Conclusion. The definition of linear phase shift networks has been generalized to 
include networks with a transfer function T(w) exp (—jkw), where T'(w) is allowed to take 
on negative as well as positive values. If the network is to be physically realizable, then 
T'(w) must be expressible by the following Fourier integral: 

ak 
T(w) = | H(u) cos wu du. 


One means of interpreting this result is that T(w) can be specified at points which are 
no closer than 2/k radians per second apart. 
A technique is presented for approximating any arbitrarily specified amplitude 
function 7;(w) by a physically realizable amplitude function 7,(w) which is given by 
k 
Tle) = | Hy(u) cos wu du, 
Jo 
where 


H,(u) = 2 | T ,(w) cos wu du “us & 
T Jo 











1960) REALIZABILITY OF LINEAR PHASE SHIFT NETWORKS 35 


If 7,(w) > 0 and is bandwidth limited, so that T,(w) = 0 if w > w, , then the unit 
impulse response of the approximating network will be positive for all values of time, 
if k < x/(2w.). If, in addition to the requirements just imposed upon 7';(w), it is also 
a monotonically decreasing function of w then the unit impulse response of the approxi- 
mating network will be non-negative if k < w/w, . Finally, it is shown that if 


[ T,(w) coswidw > 0 forall u< M, 


then, provided that k < M, any physically realizable approximate of 7',;(w) will have a 
non-negative unit impulse response. If the value of k is increased without altering 7',(w) 
then the transient response of the network will be delayed in time but its shape will not 
be altered. Thus, the value of k discussed in this paragraph can be considered to be an 
upper bound on the upper limit of integration in Eq. (7) and a lower bound on the 
magnitude of the phase slope. 
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BOOK REVIEWS 


(Continued from p. 30) 


Handbook of supersonic aerodynamics. Section 7—Three-dimensional airfoils. Produced 
and edited by the Aerodynamics Handbook Staff of the Johns Hopkins University 
Applied Physics Laboratory, Maryland, 1957. For sale by the Supt. of Documents, 
U. 8. Government Printing Office, Washington 25, D. C. 86 pp. $1.50. 


Section 7 of series first reviewed here in Vol. XVII, No. 3, p. 298. Results on those parts of lift, 
drag and pitching moment of thin wings at very small incidence which can be obtained from the Linear- 
ised Theory of steady potential flow are collected from many sources and displayed in the form of 
graphs. 

R. E. MEYER 


Applications of finite groups. By J. 8S. Lomont. Academic Press, New York, London, 
1959. xi + 346 pp. $11.00. 


This book is concerned with applications of group theory to problems of physics, and with the 
parts of group theory itself (chiefly representation theory) which are needed in these applications. 
Because of the style and the arrangement of the material, the book stands somewhere midway between 
a textbook or treatise in the usual sense and a handbook, in which a great number of definitions and 
theorems are presented without much in the way of motivation, and with almost nothing in the way 
of proof, Seemingly, the usefulness of the book will be limited, for the most part, to persons who already 
have a good general knowledge of the field, and who are in search of technical information about par- 
ticular details. 

The applications discussed include applications to crystallography, thermodynamics, wave-guide 
theory, and to nuclear, atomic, and molecular theories. As is suggested above, most of the remaining 
material centers around the representation problem for finite groups. This part of the discussion goes 
quite deep, and includes much that is certainly not widely known. There is a general bibliography 
containing over two hundred entries; and the book terminates with three appendices, one of which 
gives a twenty-five page summary of the properties of the Lorentz groups. 

It would be easy for numerous errors to lurk undetected in the dense mass of material given her 
It is only 
none of any consequence. On the other hand, an occasional poorly thought-out sentence suggests that 


proper to state, however, that the reviewer discovered only very few errors of any kind, and 






the book was perhaps written in excessive haste. If the author had reflected more deeply on the question 


i 
of what his book was intended to accomplish, and on the question of why group theory affords useful 
tools not found in classical analysis, he might have presented the subject in a more readable form, and 
appealed to a larger public. It would have been of great value, also, if he had given more specific refer- 
ences to places where proofs of the less familiar results can be found 
L. A. MacCout 
Numerical methods for nuclear reactor calculations. By G. I. Marchuk. Translated from 


Russian. Consultants Bureau, Inc., New York, and Chapman & Hall, Ltd., London, 
1959. 295 pp. $60.00. 


This is the first book to discuss comprehensively and profoundly the many interesting problems 
for the numerical analyst arising in the practical calculation of nuclear chainreactor systems. This is 
done to the exclusion of the physical background, so that questions of the validity of the various approxi- 
mations are neglected, but the author disclaims any intention of providing a comprehensive treatise 
on nuclear reactor theory. Instead, he derives the Boltzmann equation for neutron transport, presents 
approximations useful for its simplification, and deduces the finite difference equations corresponding 


(Continued on p. 60) 











HOMOGENEOUS SOLUTIONS IN ELASTIC WAVE PROPAGATION* 


BY 
JOHN W. MILES 
University of California, Los Angeles 


Summary. Busemann’s method of conical flows is formulated for two-dimensional 
elastic wave propagation. The equations of motion are reduced to either Laplace’s 
equation in two dimensions or the wave equation in one dimension, and solutions then 
are obtained with the aid of complex variable or characteristics theory, respectively. 
Special attention is paid to that class of problems in which the hyperbolic domains (of 
the two-dimensional wave equation) are simple wave zones, in consequence of which the 
solutions may be continued into the elliptic domain (of Laplace’s equation) without 
explicitly posing the boundary conditions on the boundary separating the two domains. 
The method is applied to the diffraction of P- and SV-pulses by a perfectly weak half- 
plane. 

1. Introduction. The purposes of this paper are(a) an exposition of the use of homo- 
geneous solutions’ in problems of two-dimensional, elastic wave propagation and (b) a 
complete solution for the transient diffraction of dilatational or vertically polarized 
shear waves by a perfectly weak half-plane. 

Homogeneous solutions to the wave equation appear to have been discussed originally 
by Green [1] and were developed extensively by Bateman [2], but the most powerful 
applications were initiated only much later by Busemann [3] in his method of conical 
flows. Busemann’s essential contribution was the introduction of Chaplygin’s trans- 
formation [4] to reduce the wave equation to Laplace’s equation, thereby permitting 
the use of function theory; (we note a striking adumbration of Busemann’s work in 
Donkin’s formula [5]). Busemann’s technique has been applied extensively to supersonic 
wing problems [6] and to transient diffraction by a half-plane [7, 8] and by a wedge 
(9, 10]. We also note that a somewhat earlier reduction of the wave equation in three 
dimensions to Laplace’s equation in two dimensions was effected by Sobolev [11] through 
superposition of plane waves. Sobolev applied his technique to the impulsive loading 
of an elastic half-space, but he does not appear to have emphasized either the important 
role played by homogeneity or the potential generality. 

[Note added in press. It appears that Sobolev may have given a more general treat- 
ment in Ch. XII of the 1937 Russian translation (and extension) of P. Frank and R. von 
Mises, Differential and integral equations of mathematical physics, a reference cited by 
M. M. Fridman, Dokl. Akad. Nauk. SSSR 66, 21-24 (1949)]. 

We may anticipate a homogeneous solution to a physical problem if either the data 
of that problem contain no characteristic length or if the only characteristic length 
must be derived from a parameter to which the solution must be proportional. The former 
category, in which supersonic flow past a semi-infinite cone provides an example, requires 


*Received December 8, 1958. Presented at the Institute of Geophysics Tenth Annual Conference, 
University of California, Berkeley, Dec. 13 and 14, 1957. An expanded version of this paper appeared 
as Ramo-Wooldridge Report GM-TR-0165-00350. 

1A function f(z; , 2, *** , tn) is said to be homogeneous of degree m if it may be expressed in the 


form x7" f(@2/21 » *** 5 Xn/%1). 
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no assumption of linearity, but in the latter category linearization of the boundary 
conditions generally is a necessary prerequisite’ (although the differential equation 
usually need not be linearized). We also remark that solutions to the wave equation 
may be patched together along characteristics and that homogeneous solutions then 
may be used over limited regions; the essentia] requirement is invariance under a scale 
transformation of the region in question. 

The immediate advantage of homogeneity is that the order of the governing partial 
differential equation may be reduced. This advantage proves especially great in those 
two-dimensional problems of elastic wave propagation that permit the displacement 
potentials for vertically polarized motion or the displacement for horizontally polarized 
motion to be posed as homogeneous functions of degree zero, for then we may apply 
Chaplygin’s transformation directly (as in Sec. 3 below). To be sure, we could apply 
the same transformation if the vertically polarized displacements were homogeneous 
and of degree zero, but then we should find it necessary to impose appropriate com- 
patibility relations (as in the typical conical flow problems, where the velocity com- 
ponents are of degree zero and the velocity potential of degree one; see reference [6)). 

We consider in Sec. 2 the equations of motion for homogeneous wave functions of 
degree zero, but we remark that—insofar as homogeneity of any degree may be assumed— 
the restriction to degree zero is not essential; more general results may be constructed 
via Duhamel superposition by virtue of the assumed linearity. We then go on, in Sec. 3, 
to reduce the wave equations to either Laplace’s equation in two dimensions or the wave 
equation in one dimension and to construct general solutions with the aid of analytic 
function or characteristics theory, respectively. 

Having the general formulation, we consider as an example the problem of transient 
diffraction by a half-plane when both dilatational and shear waves are present. This 
problem also has been solved by de Hoop [12], who formulated it as a Wiener-Hopf 
integral equation.” A closely related problem, the sudden opening of a semi-infinite 
crack in a previously uniform tension field, was solved by Maue [13], whose work appeals 
to the homogeneity properties discussed herein but culminates in an integral equation 
that he solves by the Wiener-Hopf technique. Maue’s results could be applied directly 
to the normal incidence on a half-plane of two, symmetrically disposed compression 
waves characterized by a step discontinuity in normal stress. 

[Note added in press. An incomplete solution based on Sobolev’s method (l.c.a.), assum- 
ing an incident wave proportional to the time-integral of (4.1) below, has been given by 
M. M. Fridman, Dokl. Akad. Nauk. SSSR 66, 21-24 (1949). He attempts only to 
reduce the problem to quadratures, the completion of which would appear to be rather 
difficult}. 

Other problems to which the present method is applicable include the action of an 
impulsive line-load on a semi-infinite solid, for which many previous solutions exist [1-4], 
and the as yet unsolved problem of P- or SV-diffraction by a wedge-shaped cavity* (the 
corresponding SH problem is essentially that solved in references [9] and [10]. 





2Problems in which the solution is nonlinearly proportional to some input parameter also may be 
possible. 

*Dr. de Hoop worked on the problem at the Institute of Geophysics, University of California, 
Los Angeles in 1956-7, His subsequent completion of the solution and the present work then were carried 
out independently. 

‘This problem is being attacked by one of Professor Knopoff’s students. 
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Equations oF MoTIon 

2.1 The displacements. We consider those (essentially two-dimensional) wave mo- 

tions in an isotropic, elastic medium for which the displacement vector u may be derived 
according to 

u= V¢4+ V x (ky) + kw, (2.1) 

where @ denotes the scalar potential for irrotational displacements (dilatational or 

P-waves) in a plane (the z, y- or r, 6-plane) normal to the unit vector k (along the z-axis), 

ky the vector potential for solenoidal displacements in this same plane (vertically 


polarized shear or SV-waves), and w a transverse displacement (horizontally polarized 
shear or SH-waves)’. It is known that ¢, ¥, and w satisfy the wave equations [15] 


oV'é = $b: ; eV t= vu, CV w= Wi , (2.2a, b, ©) 
where c, and c, denote the velocities of dilatational and shear waves, respectively, and 
are given by 


ci = (A + 2u)/p and cs = u/p (2.3a, b) 


in terms of Lamé’s constants \ and u and the density p. We also introduce the speed 
ratio y and the critical angle 0, according to 


y = cos 0. = ¢2/e, = [u/(A + 2y)]'”. (2.4) 


We now assume the boundary conditions to be such that ¢, ¥, and w must be homo- 
geneous functions of degree zero in the polar coordinates r and @ and the time ¢. We 
choose as our dimensionless, homogeneous coordinates the angle 6 and either 


€=c¢t/r or » =¢,t/r, (2.5a, b) 
in terms of which 
@ = of, 6), y = ¥(n, 6), w = w(n, 8). (2.6a, b, ce) 
The radial (uv) and tangential (v) components of displacement then are found to be 
[note that £(0/d&) = n(0/dn) = — r(d/dr)| 
u=r'(—td: + vo) and v =r" + ny,), (2.7a, b) 


where subscripts denote differentiation (as they do everywhere in this paper except in 
the following section, and, throughout, on the stress components 7;,;). 
2.2. The stresses. The Cartesian stress tensor is given by 


du; , du; 
tT; = 6; Viuct (2 + au), (2.8) 


where u; denotes a Cartesian component of the displacement with 7 = 1, 2, 3, 6;; denotes 





’We could achieve a more formal symmetry by deriving w from a second vector potential, but this 


would be unnecessarily circuitous. 
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the Kronecker delta, and the usual summation convention is implied. Transforming 
to cyclindrical polar coordinates and introducing £ and 7» from (2.5a, b), we obtain 


— oe oe (2.9a) 
Too = ur {[(y* — 2E )obcle + 2nY).0}, (2.9b) 
Tas = “Gee (2.9¢) 
To = wr *{—2Ed)eo + [(1 — 2n*)ya]e}, (2.9d) 
T,, = —pr 'tw,, and rt), = ur'w. (2.9e, f) 


We shall be interested especially in r,s, and r,, , which we also may express in the rather 
more symmetric forms 


oe 2 7 
Too = br Sl — 2n)¢, + 2nyo (2.10a) 
and 
-2 9 2 : 2 
To = ur an [—2nde + (1 — 27)y,]. (2.10b) 


2.3. The boundary conditions. We shall consider only boundaries that are either 
completely free (weak boundary) or completely constrained (rigid boundary). The corre- 
sponding boundary conditions at any point, where the components of the normal are 


mn, , are either 
7;;n; = 0 (weak) (2.11a) 
or 
u; = 0 (rigid) (2.11b) 


provided that the curvature of the boundary is finite. The boundary conditions at a 
sharp edge may be inferred from the requirement that the strain energy in the neighbor- 
hood of the edge remain bounded, which implies 


7,4; = Or”), r—0 (2.12a) 
or, equivalently, 
u, = O7'”), r—0. (2.12b) 


TRANSFORMATION OF THE WAVE EQUATION 


3.1. Wave equation in homogeneous variables. The end result of expressing V’’¢ in 
r and @ in (2.2a) and then posing the homogeneous solution (2.6a) is 


(6° — Idee + Eb: + doo = O. (3.1) 


We may obtain the corresponding equations for y and w simply by replacing & by 7. 
We remark that (3.1) is elliptic in £ > 1 and hyperbolic in ¢ < 1. If we regard &* 

and @ as polar coordinates (thereby referring all lengths to ¢,t), the elliptic and hyperbolic 

domains correspond to the interior (¢ > 1) and exterior (¢ < 1) of the unit circle (¢ = 1). 
The circle = 1 or r = ¢,t evidently represents a singular wave front, across which 
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we may expect discontinuities in the &derivatives of sufficiently high order and in the 
neighborhood of which the solutions to (3.1) cannot be uniformly valid. Fortunately, this 
neighborhood proves to be extremely small for elastic solids, and we shall rest content 
with the statement that the solutions to (3.1) must be continuous across = 1, while their 
first &derivatives may be infinite like ( — 1)" as & + 1 +. We also remark that the 
second derivatives of these solutions may be discontinuous across the characteristics 
of (3.1) in its hyperbolic domain. 

3.2. The Chaplygin transformation. We may reduce (3.1) to Laplace’s equation 


du +o. = 0 (3.2a) 
for points inside the unit circle through Chaplygin’s transformation 
s = —cosh™'t = log [§ — ( — 1)'”]. (3.3a) 
Similarly, we obtain the one-dimensional wave equation 
gee — deo = O, (3.2b) 
with 
o = cos &. (3.3b) 
It follows that we may pose the solution for ¢ in the complementary forms 
@ = RIF(a), §é>1 (3.4a) 
with 
a = 6+ cosh” &, (3.5a) 
where F(a) is an analytic function of the complex variable a; or 
¢ = F,(@@,) + F_-(@), (3.4b) 
with 
a. = 6+ cos é, (3.5b) 


where F. are arbitrary functions of the characteristic variables a. . 


-| 
i cosh € 


(€>1) 











Fra. 1. The a-plane. 


We may interpret the solution (3.4a) in terms of the mapping of the interior of the 
unit circle into the interior of a semi-infinite rectangle in the a-plane, as shown in Fig. 
1 (we distinguish the planes @ = 0 and @ = 2r in anticipation of the half-plane diffraction 
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problem to be treated in Sec. 4; if a finite sector of the unit circle were excluded by the 
physical boundaries, as in the problem of diffraction by a wedge, we could subject a to 
further transformations). We also find it convenient—especially in discussing singu- 
larities—to introduce the transformation 


z= cosa = £ cos 6 — i(? — 1)” sin 0. (3.6) 
The interior of the unit circle then maps on the exterior of the cut from — 1 to + 1, as 


shown in Fig. 2. 


~iVé2-1 sin @ 


€=1+i0 6=29 
ee . € cos @ 
-| l€=1-i0 4), 820 











Fia. 2. The z-plane. 


We may construct the characteristics of (3.5b) by drawing tangents to the unit 
circle in the clockwise (a.) and counterclockwise (a_) directions, as shown in Fig. 3. 
Any point in the hyperbolic domain (¢ < 1) then may be 





(a, ,a_) 


Fia. 3. Characteristics generated from the unit circle. 


located by the intersection of two characteristics of opposite (+) families, and we may 
designate the coordinates of this point as (a, , a_); if a, = a_ the point is on the unit 
circle (¢ = 1), while if a, = a. + = the point is at infinity (¢ = 0). 
3.3. Application to shear waves. We may apply the results of Secs. 3.1 and 3.2 
equally to ¥ and w if we replace £ by 7; thus, 
y = R 1G(8) (3.7) 
with 
B= 6+ icosh” 7, n> 1 (3.8a) 
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or 
B. = 6+ cos 7, 9 <1, (3.8b) 


and similarly for w. 

3.4. Simple wave zones. A simple wave zone is a region in which the solution to a 
hyperbolic equation depends on only one of the characteristic variables. The simple 
wave zones that we encounter in the present study are bounded by a characteristic, along 
which they are adjacent to a region of zero disturbance’; by an arc of the unit circle, 
along which they are adjacent to an elliptic region; and by a boundary on which appropri- 
ate conditions are prescribed, as shown in Fig. 4. 

We infer from the foregoing definition and from the required continuity of ¢, = 
Ri F’(a) that the solutions given by (3.4a) may be continued into a simple wave zone 
simply by continuing @ as the variable (either a, or a_) that is constant along the 
characteristic boundary; in so doing, we may drop the + or — subscript. We also may 
continue implicitly the prescribed boundary condition along the single family of character- 
istics to the unit circle and thereby seek the solutions in the elliptic and hyperbolic 


-lp_ 
8+cos €=a, 


6+cos '€= a<o, 


ay 





Fic. 4. Simple wave zone bounded by an a, characteristic, an arc of the unit circle (0 < 6 < a4), and 
a line (@ = 0) of specified boundary condition. 


domains of (3.1) simultaneously (whereas in a more general configuration it would be 
necessary to obtain the solution first in the hyperbolic domains in order to specify 
boundary conditions completely around the boundary of the elliptic domain). 

We consider, as a particular example, a dilatational disturbance originating at r = 0 
and ¢ = 0 and moving out along a free boundary @ = 0. In order to satisfy the boundary 
conditions there, shear waves generally must be excited by the dilatational wave front; 
in accordance with Huygens’ principle, each of these may be regarded as a cylindrical 
wave having its center at r = c,r and 6 = O and having a radius c, (¢ — 1), 
where 0 < + < ¢. The envelope of these waves then consists of that portion of 
the characteristic 6 + cos”' » = 6, intercepted between the wave front r = c.f and the 
free boundary, as shown in Fig. 5. Such envelopes (which also arise in problems, elastic 


Co(t-t) 





a 


Fic. 5. Simple wave zone for shear waves excited by dilatational wave front moving along @ = 0. 











c,t 





‘It follows from characteristics theory that a simple wave zone is the only type of hyperbolic region 
that may adjoin a region of zero (or uniform) disturbance. 
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or acoustic, where two media of different wave speeds are contiguous) sometimes have 
been designated as head waves; however, following the terminology of characteristics 
theory in supersonic flow [16], we shall designate them as Mach lines or Mach envelopes. 

3.5. Transformation of the displacements. We consider only the radial and 
tangential displacements u and v; the transverse displacement w is given directly by a 
solution of the type (3.7) and requires no further transformation. Substituting (3.4a) 
and (3.7) in (2.7a, b), we obtain 


u = r'RI[—i( — 1)'? F(a) + G(8)] (3.9a) 


and 
vy =r RI[F’(a) + in(n’ — 1)7°'”G"(8)] (3.9b) 


for points in the elliptic domains; if a is continued as a. we have only to replace 
i(e? — 1)°°”? by + (1 — &)°™ and similarly for 8. 

The results (3.9a, b) may be continued into simple wave zones, but we shall find it 
more convenient to utilize the fact that the P- and SH-wave displacements must be 
normal and parallel to their respective characteristics (see Fig. 6). Designating these 
displacements as u* and v® (v® is measured positive in the direction of increasing @ for 
both + characteristics), we find 


I 


us = +t? —ct)'"o., a= O+cos'é (3.10a) 


and 
Y= Fr — Bt’) y, , 8 = 0+ cos 7. (3.10b) 


We emphasize that the total displacements generally contain both P- and SH-com- 
ponents and that the hyperbolic domains of these components do not coincide. 

3.6. Solution procedure. The foregoing developments suggest the following, general 
procedure for the solution of those problems that admit homogeneous functions of 
degree zero and in which the hyperbolic domains may be identified as simple wave zones. 

a. Pose the solutions in the form of (3.4a) and (3.7). We shall find it expedient, in 
so doing, to separate out prescribed components, such as incident waves [e.g., (4.2a, b)]. 

b. Convert the initial conditions to conditions on the wave fronts [e.g., (4.4a, b)]. 

c. Calculate the displacements and/or stresses in forms similar to (3.9a, b) and 
impose the prescribed boundary conditions [e.g., (4.6a, b)]. 

d. Express the coordinates that vary on the boundaries (¢ and/or 7 in Sec. 4, where 
the boundaries are prescribed by fixed @) in terms of a and £ there (e.g., § = cos a on 
§ = 0) and substitute in the boundary conditions to obtain equations that contain only a 
and/or B. 

e. Express the prescribed functions (representing applied loads, displacements, or 
incident waves) in terms of a or 6 on the boundaries. We shall find (Appendix A) that 
the resulting expressions are characterized essentially by their singularities (generally 
poles) but contain arbitrary functions. 

f. Determine the relation between a and 6 on the boundaries (e.g., cos 8 = y cos a 
on 6 = 0) and eliminate (at least implicitly) either a or 6 from the boundary conditions. 
We remark that either the same or different relations between a and 8 may be obtained 
on different boundaries; thus, we shall find only a single relation in Sec. 4, but two 
different relations would be obtained in the problems of P- and SV-wave diffraction 


by a wedge. 
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ON 








Fia. 6. The velocities in the hyperbolic domains. 


g. Determine functions that satisfy the boundary conditions. We find that, in the 
problems to be considered subsequently, only the first derivatives of F and G enter the 
boundary conditions; accordingly, we may continue these conditions into either the a or 
8 domains (see f above) and solve for F’ and G’. If the form of the boundary conditions 
is different on different boundaries the resulting expressions for F’ and G’ must be de- 
veloped as superpositions of the algebraic solution of the different conditions. 

h. The solutions so obtained will contain the arbitrary functions introduced in e, 
and we must determine them in such a way that the final solutions satisfy the initial 
conditions, exhibit singularities on the physical boundaries only as physically appropriate 
(e.g., stress may be infinite like r~'”* at an edge), and are completely analytic in the 
interior of their elliptic domains. We find that the first two of these conditions may be 
met in a relatively straightforward manner (indeed, the initial conditions are likely to 
be satisfied automatically once the wave front geometry is determined). We also find 
that solutions satisfying the third condition of analyticity often may be obtained 
by inspection, but in See. 4 the removal of spurious singularities introduces appreciably 
greater complexity into the analysis. We find that the undesired singularities may be 


factored out through Cauchy integral-theorem representations, as in the solution of 


Wiener-Hopf integral equations. 


DIFFRACTION BY HALF-PLANE 
4.1. Formulation. We consider (see Fig. 7) the P-wave 


¢* = Hit — cos (0 — 4,)] (4.1) 
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to be incident on the weak half-plane 6 = 0, 27; H denotes Heaviside’s step function 
[H(x) = 0,1 asx <, > O] and 6, the angle of incidence. It would be equally simple, in 
principle, to provide for the simultaneous incidence of an SV-wave; this would almost 
double the length of the subsequent equations, however, and we shall rest content with 
sketching in the results by analogy (see Sec. 4.5). 








Fic. 7. Diffraction of P-wave by half-plane. 


Our problem evidently admits no characteristic length other than that defined by 
the amplitude of the incident wave (which we take to be unity); accordingly, we may 


pose solutions in the form 


¢=¢ +RIF(a) and y = RIG(8). (4.2a, b) 

We remark that F(a) and G(8) will include specularly reflected P- and SV-waves emerg- 
ing at the angles 27 — @, and 2x — 6, , respectively, where 

cos 6. = y cos 6, . (4.3a) 


If 0 < 6, < x/ 2¢' constitutes the entire disturbance for ¢ < 0, but if r/2 <0, <7 


these specular reflections, as given by (4.24b, c) below, must be included for ¢ < 0 as 
well as ¢ > 0. We also remark that (4.3a) may be generalized to read 


cos 8B = y cosa, 6= 0, 2r, (4.3b) 


since 7 = yé on the half-plane; (4.3b) then yields (4.3a) at the point of intersection of 
the plane-wave fronts (see Fig. 7). 

Turning to the scattered waves, we first note that, since the P-wave disturbance from 
the edge at r = 0 originates at ¢ = 0 and travels outward with speed c, , the region of 
influence of this edge consists of the interior of the circle £ = 1, which we then may 
identify as the elliptic domain for F(a). Similarly, the elliptic domain for G(8) is the 
interior of the circle 7 = 1; but, in accordance with the discussion of Sec. 3.4, the SV- 
wave region of non-specular scattering also comprises the simple wave zones bounded 
by the Mach envelopes @ + cos‘ 7 = 6, and 6 — cos 'n = 27 — 6, . Summing up, we 
designate the various zones according to (see Fig. 7): 


I: incident wave zone, 
II, : P-wave zone of specular reflection, 
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II, : SV-wave zone of specular reflection, 
III: shadow zone, 

IV, : P-wave scattering zone (~ > 1), 
IV, : SV-wave scattering zone; 


this last zone comprises both elliptic and hyperbolic domains. 

The initial conditions dictate: (a) ¢ = ¢° on that portion of the P-wave scattering 
circle £ = 1 intercepted between the incident wave front and the specularly reflected 
P-wave front and (b) ¥ = 0 on the Mach envelope 6 + cos‘ 7 = 6, and on that portion 
of the SV-wave scattering circle intercepted between this Mach envelope and the 
specularly reflected SV-wave front (note that 6. > 6. for all 6,). Imposing these con- 
ditions on (4.2a, b), we obtain 

Ri F(6) = 0, $< 0<2r — O, (4.4a) 
and 
RI G(6) = 0, 6.5 0<2r — 6, (4.4b) 


where the end-point 6 = @. comprises the aforementioned Mach envelope. 
The boundary conditions require 7,, and 79, to vanish on the half-plane; integrating 


(2.10a, b) from » = 0 (i.e., either = Oorr = @) yields . 
(1 — 29')o, + 2ny. = 0 (4.5a) 
and 
—2nds + (1 —2n)¥,=0, O=0, QO. (4.5b) 
Substituting (4.2a, b) in (4.5a, b), we obtain 
y'(1 — 2n°){6(E — cos 0) + RIfi(e? — 1)7'?F’(@)]} + 2nRIG’(8)] = 0 (4.6a) 
and 
—2n{—sin 6,5(E — cos 6,) + RI[F’(a)]} + (1 — 2m’) (4.6b) 
-RIfi(n? — 1)7'?G’(8)] = 0, 6 = 0, 2r. 
4.2. Solution for F’ and G’. We now rewrite (4.6a,b) in terms of a and # and repre- 
sent the delta functions as in Appendix A to obtain 
Rifese a cos 28F'(a) + 2y cos BG’(8)] = RI[A(cos a)/ix(cos a — cos 6,)], (4.7a) 
Ri{—2 cos BF’(a) + ese 8 cos 28G’(8)] = —RI[B(cos 8)/ix(cosa — cos 6,)], — (4.7b) 


A(cos a) = +cos 20, , a= 4 , (4.8a) 
2Qr _— 6; 


so ™ (4.8b) 


2r — 6, 


B(cos 8) = +ysin 20, , 


and 


Im A(é) = Im B(y) = 0. (4.8¢c) 
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Our choice of arguments in A and B is dictated merely by convenience; in this connec- 
tion, we note that 8 = @, is equivalent to a = @, [see (4.3a, b)] in (4.7a, b). 

We may continue (4.7a, b) into either the a- or 6-domain by eliminating either 
6 or a, respectively, through (4.3b) and removing the R/ operators [i.e., the equations 
so obtained yield solutions that satisfy the boundary conditions (4.3b) and (4.7a, b)]. 
Solving the resulting equations, we obtain 


F’(a) = sin a(A cos 28 + yB sin 28)/ix(cos a — cos 6,)D(y cos a 4.9a) 
and 
G’(8) = sin B(yA sin 2a — B cos 28)/imr(cos a — cos 6;)D(cos 8), 4.9b) 
where 
D = (2y’ cos? a — 1)? + 47’ sin a cos’ a(1 — 7’ cos’ a)! 4.9¢e) 
= cos 28 + 4sin 6 cos’ B(y’ — cos B)'”’, 1 .9d) 


with cos 8 and cos a defined by (4.3b) in (4.9a) and (4.9b), respectively. 

The denominator D(y cos a), as given by (4.9c), has branch points at cos a = -+ 1/y 
and (at least for normal values of y) zeros at cos a = + £, , where 

ER = C, Cr; $.10a) 

and cp denotes the wave speed for Rayleigh surface waves. Similarly, D(cos 8) has 
branch points at 8 = 6, and  — 6, (cos 8 = + cos 6.) and zeros at cos 8B = + ne, 
where 
‘4.10b) 


Nr = Yor = Co/Cr- 


tr 


We choose the branch cuts so that they do not enter the domains of regularity, — > 1 
and 7 > 1, and define the radicals to be positive and real for cos a = 0 and cos 6 = 0. 
We find it convenient, in the subsequent determination of A and B, to introduce the 


notation 
Zz, = cosa and z, = cos®, (4.11la, b) 


@1 
and to examine the singularities of F and G in z,- and z,-planes. The elliptic domains 


map onto the z,- and z.-planes cut from — 1 to + © (the cuts from — 1 to 1 correspond 


to £ = 1 and n = 1, while the cuts from + 1 to + © correspond to the half-plane separat- 


VV 


ing 6 = 0 and 6 = 27), and we require 








dF _ _ Fa) _ Aa)(2y’zi — 1) + 27° Boa)a(l — ¥'2)'” 112 
dz, sin a (—in)(z; — cos 6,)D(y2;) oo 
and 
dG _ _G"(B) _ 2A(z./y)ze(y* — 22)? — yBle.)(222 — 1) {.12b) 
dz, sin 6 (—1t71)(z. — cos 62) D(z-) a 


to be regular in these cut planes. We also require 


dF /dz, = O(z;'”’), Zz, © (4.12¢) 
and 
dG/dz, = O(z;"”’), 22 —> © (4.12d) 
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in order to meet the edge singularity requirements of (2.12), which imply (see Appendix 
B) F’(a), G’(8) = O(r-”*), r 0. 
The requirement that dF /dz, be regular in the z,-plane cut from — 1 to + ~ eee 


that A(z,) and B(yz,) must cancel the singularities exhibited by D(yz,) at z, = — 1/y 
and z, = — &p [see (4.9c) and discussion following]; if we satisfy this requirement, it 
will follow that A(z./y) and B(z2) will cancel the singularities of D(z,) at z2 = — 1 and 
Z. = — ne, although only the latter cancellation may be posed as an a priori require- 
ment (indeed, dG/dz, still will possess a branch point at z. = — 1). We can achieve 
these cancellations by factoring D according to 

D(z) = D(z) D_(2), (4.13) 


such that D.(z) have singularities only in Riz 2 0. 

At this point, our analysis parallels that of both de Hoop [12] and Maue [13], for 
D(z) enters the transform kernel of the Wiener-Hopf integral equation and must be 
factored in the same way; we obtain’ 





D.(z) = [2(1 — y*)]'"(ne F 2)/L(+ 2), (4.14) 
where 
- Lf’ xe) } 
Az) = ex {1 f >a ‘ 
L(z) Pi) poet (4.15) 
and 
Ris. gryi/ts? 1/2 
x(@) = tan” k ‘Sook 5 tit Sool & |,o éxs% (4.16) 
c= F 2 


is the phase angle of D(z) for points on the bottom of the cut from ¢ = 0 tot = y. The 
path of integration in (4.15) must be indented over or under z if z tends to the real axis 
interval (y, 1) from Im z < 0 or > 0, respectively, and we then may rewrite (4.15) 
according to 


0O< 6< @ 
( . —_ ix (eos 6) aut at 
L(cos 6) = Le(cos be , or — 0. < 0 < 2n’ (4.17) 
where L, denotes L based on the Cauchy principal value of the integral in (4.15)—viz., 
‘j = aay x(z)/* { 1 x(t) _ x(x) } 
it) = xp 4- ————-— < # < i. : 
A-E= exp 4) pa hy ses (4.18) 


Numerical values of x(z) and L(z) have been calculated* for y = 37’ (A = uw) and 
are plotted in Figs. 8 and 9. We remark that these results have application to other 
problems in elastic wave propagation (e.g., Maue’s results, [13]). 

Having the result (4.13), we can cancel the undesired singularities in D by including 
the factors D_(yz,) and D_(z,) in A and B, respectively. Then, since D_(yz,) takes the 
same values at a = 6, and 2x — @, (z, = cos 6,) and is O(z,) at infinity, we require an 
additional factor in A that takes opposite values on a = 6, and 27 — @, in consequence 
of (4.8a), is O(z,~'”*) in consequence of (4.12c), and is regular in the cut 2,-plane; we 





7The details of the factorization may be found in Ref. [13]. The results were derived independently 
by the writer and therefore provide a check on those of Refs. [12] and [13]. 
8] am indebted to Professor Knopoff for aid with these and the subsequent calculations. 
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Fic. 8. (2/x)x(z), as given by (4.16) for y = 3-1/2, 


can satisfy these requirements with the factor (1 + z,)~’”*. Similar requirements apply 


to B, but B(yz,) also must remove the branch point at z, = — 1/y in the numerator of 


dF/dz, ; we can satisfy all of these requirements with the factor (1 + z,)~'””. Introducing 
the normalization factors dictated by (4.8a, b), we conclude that 


A(cos a) = cos (6,/2) cos (262) sec (a/2) D_(y cos a)/D_(cos 62) (4.19a) 
and 
B(cos 8) = vy sin (26,) cos (6,/2) sec (8/2) D_(cos 8)/D-_(cos 62). (4.19b) 
Substituting (4.19a, b) in (4.9a, b) or (4.12a, b), we obtain the final results 


{C,, sin (a/2)(2y’ cos’ a — 1) + yC2, sin (2a)[(1 — y cos a)/2]'”"} L(y cos a) 








, 7 Pon Li 
“ee = imy(—ép — cos a)(cos a — cos 8,) : 
(4.20a) 
G’(8) = {C,, sin (28) [y(y —_cos B)/2y" — C,, sin (6/2) cos (28)}L(cos B) ; (4.20b) 
it(nr — cos B)(cos B — cos 62) 
C;; = cos (6,/2) cos (26.)C,(cos 62), (4.21a) 
Cx = y’ sin (26,) cos (02/2)C,(cos 62), (4.21b) 
and 
C,(cos 62) = L(—cos 62)/(1 — y’)(nk + cos 6), (4.21c) 


where L is defined by (4.15), a by (3.5), 8 by (3.8), y by (2.4), &g and n»z by (4.10a, b), 
and 6, by (4.3a). We observe that C,, and C., are real for all 6, and that C,, = 
—(1 — y’)"’” and C., = 0 for normal incidence (6, = 2/2). We also remark that the 
second (1) subscripts on C,, and C., identify the excitation as a P-wave and that the 


results (4.20a, b) agree with those obtained by de Hoop [12]. 
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Fic. 9. The function L(x), as given by (4.15) for —1 < x < y and by (4.18) fory < z < l withy = 3-12. 


4.3. The scattered wave fronts. We consider first the plane wave zones, where the 
solutions may be obtained directly from (4.8a, b) and (4.9a, b). Outside of the P-wave 
circle £ = 1, both a and 8 are real, and the only contributions to Rl F’(a) and RI G’(8) 
must come from the poles at a = 0,,a = 27 — 6, , and 6 = 27 — @, (we remark that 
the numerator of G’ vanishes at 8 = 6, , thereby cancelling the corresponding zero in the 


denominator); these poles yield 


RI F’(a) = sin 6,4[~ cos 6 — (1 — &*)'” sin 6 — cos 8,] in IIT, (4.22a) 
RI F(a) = R,, sin 6,4[§ cos 6 + (1 — &)'” sin @ — cos O,] in II, , (4.22b) 
and 
R1G’(8) = Rx sin 6,6[n cos 6 + (1 — 7’)'” sin 6 — cos 6] in II, , (4.22c) 
where 
R,, = (—cos’ 26, + 7 sin 26, sin 20,)/D(cos 62), (4.23a) 


R., = 2y’ sin 26, cos 26,/D(cos 62), (4.23b) 
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and 

D(cos 6.) = cos’ 26, + y’ sin 20, sin 26, . (4.23¢) 
Integrating (4.22a, b, c) and modifying the argument of the resulting step functions 
(only the locus across which this argument changes sign being important), we obtain 


the scattered waves 


¢’ = RI F(a) = —H[6, — (6 + cos §)] in III, (4.24a) 

¢° = RI F(a) = R,,H[(@ — cos” &) — (2r — @,)] in II, , (4.24b) 
and 

v¥ = RIG(6) = RaH[(@ — cos' n) — (24 — 0,)] in II, . (4.24¢) 


We have, then, the anticipated results that ¢* cancels ¢* in the shadow zone (III), and 
that ¢* and y in the zones of specular reflection (II, .) agree with the known results [17] 
for reflection from an infinite, plane boundary. 

Turning to the P-wave front, £ = 1—, we find that F’(a) is imaginary there and 
yields only a radial velocity [sce (3.9a, b)] that tends to infinity like (#? — 1)~'” and 
then drops discontinuously to zero (except at 6 = @, and 2x — @,) for ¢ = 1+. Setting 
¢ = ¢,t/r in the (# — 1)~'” factor and ¢ = 1 elsewhere in the first term in (3.9a), we 
obtain 


u— (cit? — r’)7'”"F,,(8), r—ct—, (4.25a) 


where 
fiu(@) = Ri[—iF’(8)]. (4.25b) 


We may calculate the angular distribution function directly from (4.20a), since iF’(@) is 
real; the result evidently will be singular at 6, and 2x — @, in consequence of the plane 


waves of (4.22a, b). In the special case of normal incidence we obtain 


_ sin (6/2)(1 — 2y° cos 6)L(y cos @) . nal (4.26) 


2\1/2 nl 7 j @&e 
my(l — vy)’ cos (Ee — cos 0) ’ 








fi:(9) = 


We find a rather more complicated behavior at the SV-wave front (yn = 1) in conse- 
quence of the simple wave zones in r > c,t. We obtain for the (SV-contribution to the) 
tangential velocity just inside the semi-circular wave front 


v— (et — ry)” f.,(8), r=cCt—, (4.27a) 


where 

foi(0) = RI[iG’(6)], (4.27b) 
but this component vanishes on r = c,¢ + only for 0, < 6 < 2a — 0, . The remaining 
intervals are in the simple wave zones, where (3.10b) yields 


Ye = (r? _ ct’) "920 + cos” 7), (4.28a) 
where 


<= ¢ < 
Jo1( 8) = FRIG’(6)], ¢ 0. e ( 1 .28b) 
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We then may obtain v on r = ct + by setting cos~' » = 0 in (4.28a), but we emphasize 
that the resulting approximation is not uniformly valid in the neighborhoods of 8 = @, 
and 2x — 6. , where @’ vanishes like (cos 8 — cos 6,)'”* in consequence of the boundary 
condition (4.4b). We may calculate f., and gs; directly from (4.20b), making allowance 
for the fact that L (cos 8) is real only for 0. < 6 < 24 — 0, ; in the special case of normal 


incidence we obtain 
fo:(0) = —h2,(0)L(cos 8), 65 0S 27 — 6, (4.29a) 


—h,,(0)L¢(cos 6) sin [x(cos 6)], 


II 














0<6< 86, or 2r — 6. 5 0 < 2z, 6, = 2/2, (4.29b) 
fi 
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Fic. 10. The distribution of radial velocity along the P-wave circle 


(r = c, t —) resulting from normal 
incidence of a P-wave pulse; see (4.25a, b) and (4.26) with y = 3-1/2, The insert is a polar plot. 
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where 
A. i/3 10s 8 — » |!/2 
(@) = + ( ar) me cs a (4.29¢) 
; T = y (nr — COS 6) 
and 
Joi(@) = ho,(0@)Le(cos 6) cos [x(cos 6)], 6, = 2/2. (4.30) 
-| 
79=cos (1/3) 
5 / 
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4 
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0.4 T rig /3) 
8 = cos—' (I//3) 
yl 
| 
| || 
1 \821 
Ok 1 | | 
NTT? 
| ) | | 
-0.2+ a ee 
fo, | | fa) 
| 
| 
| ‘| 
-0.4 | | 
V 
| 
| 
-0.6 . 
180° 150° 120° 90° 60° 30° o° 
@ 


Fic. 11. The distribution of the tangential velocity along the S-wave circle resulting from the normal 
incidence of a P-wave pulse; see (4.27)-(4.30) with y = 3-1/2. The insert is a polar plot. 


The numerical values of f,,(0), f2:(9), and gz,(@), as given by (4.26), (4.29), and 
(4.30) with y = 37”, are plotted in Figs. 10 and 11. We emphasize that the discontinuity 
in f;, at @ = 0; = 2/2 is a direct consequence of that in ¢’. 

4.4. Surface displacements. We may simplify the expressions for the surface dis- 
placements by introducing (4.3b) in (3.9a, b) to obtain 


u =r 'RI{ cot aF’(a2) + G’(B)} (4.31a) 
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and 


v =r 'RI{F’(a) — cot BG’(8)}, 6 = 0, 2r. (4.31b) 


Substituting (4.20a, b) in (4.3la, b), we obtain 





7) [lv/2)C1, see (a/2) cos (28 — a) + C,, sin (8/2)]L(cos 8) 
latin imr(nr — cos B)(cos 8B — cos 62) (4.32a) 
and 
és | yOu sin (a/2) + (y/2)Cx sec (8 /2) cos (28 — a)]L(cos 8) (4.32b) 
( imr(nr — cos B)(cos B — COs 8.) J 


We may reduce (4.32a, b) further if r < ct, for then the first and second terms in 
the numerator of (4.32a) are real and imaginary, respectively, and conversely for (4.32b). 
The real terms in the numerators make no contributions to the displacements except 
at the surface wave pole, where they yield delta functions. Noting that 


y cos (28 — a) = cos 28/2 cos B, n = OR (4.33) 


in virtue of the evanescence of D(n,), we obtain 














— Cis(2nr — I) L(ne) 6Cxt ee _ Corl *Cot — r)'?L(cot/r) 
2°7(1 + En)'*ne(ne — COS 02) 2" nR(Cot — r Cos O2)(Crt — 1) (4.34a) 
and 
7 =-_Ca(2m = 1) Ln) deat — r) Pa __ Cyr" 2(e,t — r)'?L(cot/r) 
2°/7(1 + ne)'’’nk(ne — COS 02) ~ 42'*Eg(e2t — 7 COS 02)(Crt — 7)’ (4.34b) 
= he r < ef. 


4.5. Incident SV-wave. We now assume an incident SV-wave in place of the incident 
P-wave of (4.1). Replacing (4.1) and (4.2a, b) by 


@¢=RIF(a) and y = H[n — cos (6 — 6.)] + RI G(8), (4.35) 


imposing the boundary conditions (4.5a, b), and then proceeding as in (4.6)-(4.20), 
we obtain solutions identical with (4.20a, b) if C,, and C., therein are replaced by 


C,. = cos (0,/2) sin (20,)C, and C.. = —cos (62/2) cos (26.)C%o . (4.36a, b) 


We emphasize that, although @, and @, still are related by (4.3a), they now have 
somewhat different meanings—viz., 0. is the angle of incidence and has an ad- 
missible range (0, 7/2),” while 6, , formally defined as the angle of reflection for the 
P-wave, can be real only if 6, lies in (@, , x/2). If 6, does lie in this interval the wave-front 
geometry differs from that of Fig. 7 only in that the incident wave zone (I) and the shadow 





°No physical significance can be attached to +/2 < 62 < 7, for the incident SV-wave then would 
be preceded by a P-wave for all ¢ < 0, and non-plane, scattered waves also would exist for all ¢ < 0. 
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zone (III, , say) will be bounded by those parts of 8. = 27 — 6, and 8, = @,—rather 
than a. = 27 — 6, anda, = 6,—lying outside of the P-wave circle (£ = 1). We then find 
vy’ = RIG(8) = —H[é@. — (6+ cos” »)] in III, , 1.37) 
@ = RI F(a) = R,.H[(@ — cos’ §) — (2x — 6,)] in II, , (4.37b) 
v’ = RIG(8) = R..H[(@ — cos” n) — (2x — 8:)] in II, , 1.37¢) 
R,. = —2sin 26, cos 26,/D(cos 6.2), (4.38a) 
and 
Ra = FR, , 6. < 0. < 2/2, (4.38b) 


where R,, and D are given by (4.23a, c) for given @,,. . The results (4.37b, c) agree with 
those for reflection from an infinite, plane boundary [18]. 

If 6, lies in the interval (0, @.) the zones of specular reflection no longer exist. The 
angle 6, then is complex, cos 6, is real but greater than unity, and the pole at a = 27 — @, 
(F’ remains regular at a = @,) lies in the P-wave circle on the illuminated side (@ = 27) 
of the half-plane, where it represents a surface wave moving with the dilatational wave 
speed c, . The poles at 8 = 6, and 2x — @, lie inside the simple wave zones bounded by 
6 = 0,and 8 = 2nr — 86, , repectively, and represent singular components of the complete 
disturbance there, although we remark that the pole at 8 = 6, just cancels the incident 


wave in the angular interval (0, 6). 
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Fig. 12. The distribution of radial velocity on the P-wave circle (r 
incidence of an SV-wave pulse; see (4.39) with y = 371/83, 





= ¢, t—) resulting from normal 
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Turning to the non-plane wave fronts, the results (4.25), (4.27), and (4.28) remain 
valid if F’ and G’ are calculated on the basis of (4.36a, b); we may distinguish the results 
by writing fiz , for , and go2 in place of fi: , for , and go, . In the special case of normal 


incidence we obtain 





















































1 2 \'” sin 6(1 — 7 cos 0)'”L(y cos 6) rT 
( =o ( :) : = <= . 
10) = -2(7=a pee 12=5, (4.39) 
= ‘= — 9 4 < < —= 
foo( 8) hoo(@)L(cos 8), 6.< 0S 24 — 0., (4.40a) 
= —h,(6)Lc(cos 8) cos [x(cos 4)], 
0<6< 06. or 2x — 6, < 0 < 2x, 6, = 2/2, (4.40b) 
where 
be a 1 sin (0/2) cos 26 
hoo( 8) = x(1 Ae yy cos (nr — cos 6) ’ (4.40c) 
lte2| 
8 
ty 
6 : 
@=cos' (I//3) 
! 
| 
4 
| 
foo | 
2 
| 
i 
0 
\ yy, 
foo / 
922 
-2 
-4 
= ° ° ° ° ° 
180 (35 30 45 0 
8 
Fic. 13. The distribution of tangential velocity on the S-wave circle resulting from normal incidence 


of an SV-wave pulse; see (4.40) and (4.41) with y = 371”, 
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and 
Jo2(0) = —h22(6)Lc¢(cos 8) sin [x(cos @)], 6, = 2/2. (4.41) 


The numerical values of f,2(@), fo2(@), and g22(6) are plotted in Figs. 12 and 13. 

4.6. Rigid half-plane. Assuming an incident P-wave and solutions in the form of 
(4.2a, b), we find that the requirement that u and v, as given by (2.7a, b) and (3.9a, b), 
vanish on the half-plane yields [see (4.7a, b)] 


Ri{ cot aF’(a) + G’(8)} = cos 6,6(cos a — cos 9,), 1.42a) 
and 
RIU{F’(a) — cot BG’(8)} = sin 6,6(cos a — cos 4), 6 = 0, 2r. (4.42b) 


These equations may be solved as in (4.7)—(4.20), but the results are of rather limited 
interest, and we note here only that the counterpart of D is 


D=ycos(a — —) = 2+ (1 — 3)’ -— 2)”, (4.43) 


which has no zeros in the cut plane (s). 


APPENDIX A 
Representation of 5(x). We require f(z), 2 = x + ty, such that 
Ri f(z) = 62), y = O+. (A.1) 


Noting that 


“f.5 Sit 1 
6(x) = +- lim (4), (A.2) 
T 0+ 2 + Y 


we find that (A.1) is satisfied by 


f(z) = a(z)/inz, (A.3a) 
where a(z) has no other poles on the real axis and satisfies 
Im a(x) = 0, (A .3b) 
and 
aja +120) = ¥1. (A.3¢e) 


We remark that this result also could have been deduced from Cauchy’s residue theorem, 
which implies that f(z) must have a simple pole with residue + (i)~' at the origin and 
be real everywhere else on the real axis. 


APPENDIX B 
Asymptotic behavior at edge. We require the asymptotic behavior of the displace- 
ments given by (3.9a, b) in conjunction with the solutions of (4.9a, b). Expanding the 
radicals in (3.9a, b), we obtain 


u = r ORI —i[F'(a) + iG’(8)] —HE?F'(a) + OF 'F’)} (B.1a) 
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and 
v =r 'RI{[F’(@) + iG’(8)] + 4in?G"(8) + O(n “*G’)}. (B.1b) 
Now, 
cos a ~ te~**[1 + O(F*)] (B.2a) 
and 


cos 8B ~ ne~*[1 + O(n-”)], (B.2b) 


so that cos a and cos 6 are asymptotically related as in (4.3b). Substituting (B2a, b) in 
(4.9a, b), we find that 


F’(a) + 1G'(8) = O€ °F’, 1°G’). (B.3) 
It then follows from (B.1a, b) that 
u,v =r 'O(¢ °F’, 1 °G’) = OCF’, rG’). (B.4) 


Imposing the requirement (2.12), we obtain 
F',@" = 09") = 0G", 1°). B.5) 
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BOOK REVIEWS 


(Continued from p. 36) 
to the differential or integral equations of the problem. Amongst others, the diffusion and the diffusion- 
age approximations are discussed in detail. The difference equations are then carefully discussed and 
solved, and there is a chapter on pertubation theory. 

The book should be of great interest not only to specialists in nuclear reactor calculations but to 
anyone interested in the numerical solution of partial differential equations. Unfortunately, the price 
of $60.00 will take the book outside the range of most. 

WALTER FREIBERGER 


Elementary matrix algebra. By Franz E. Hohn. The Macmillan Co., New York, 1958. 
xi + 305 pp. $7.50. 
This book fulfills a long-felt need for an introductory textbook on modern matrix algebra which 
is suitable for students in all fields of applied mathematics. In particular, it provides admirably the theo- 
*kground for a course in the numerical methods of linear algebra. It is exceptionally well 


retical bax 
ith utmost clarity. There is a wealth of worked and unworked examples. The 


organized and written w 
contents cover: rules of matrix algebra; determinants; notions of inverse, rank, equivalence; linear 
dependence; aces; linear, unitary and orthogonal transformations; characteristic equation 


of a matrix; b 





*, quadratic and hermitian forms. 
WALTER FREIBERGER 


Mathematics for communication engineers. By 8. J. Cotton. The Macmillan Co., New 


York, 1959. x + 245 pp. $7.50. 
In the preface it is stated that the book is written to bridge ‘‘the gap between the level of mathe- 


matics reached on many (undergraduate) courses,..., and that required to read intelligently many 


of the articles in professional journals or to tackle the many problems facing engineers-particular]h 
communications engineers.’”’ There is a definite need for such a book. The reviewer believes this book 
will be quite useful to practicing electrical engineers needing a lucid and fairly precise discussion of 


| 
much of the mathematics usually studied by seniors or first year graduate students in electrical engi- 
neering. Because of the wide range of subjects considered, and its high level of presentation, the book 


would probably be valuable as a reference for other readers. 





A summary of the topics included may give an idea of the scope of the book; review ulgebra 
and differential calculus, determinants, Maclaurin and Taylor series including a short discussion oi 


inctions, integral calculus with schemes for evalua 


convergency and approximation, trigonometric { 
of linear 


certain types of integrals, a summary of differential equations emphasizing the solutior 
equations, a brief treatment of transmission lines and filters, theory of functions o complex 
variable including Cauchy’s integral formula, contour integration and conformal mapping, Fourier 
series and the Fourier integral, the Laplace transform and Heaviside’s expansion formula, Vector 
Analysis and the derivation of Maxwell’s equation and, in a problem, the derivation of the wave equation 
for perfect dielectrics, a brief mention of Legendre functions, Bessel functions and the Gamma function 
including the useful expansion of sin (z sin @) in an infinite series of Bessel functions, and finally, a 
chapter on probability and statistics containing Poisson’s approximation to the binomial law but 
omitting the normal distribution. The author has not attempted a complete discussion of each of the 
above subjects but has managed to include most of the areas which have been found useful up to now 
in electrical engineering. Because of this selection, the book may not be entirely satisfactory as a text- 
book for a course in mathematics. The mathematics presented is at least as sound as that given in the 
better, recent engineering or engineering mathematics textbooks; probably as important to the practicing 
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AMPLITUDES OF OSCILLATIONS GOVERNED BY A MODIFIED 
VAN DER POL EQUATION’ 


BY 
K. KLOTTER? anp E. KREYSZIG: 


1. Introduction; systems treated. The classical example of a system exhibiting 
self-sustained oscillations is the one described by the differential equation 


q’ — éq (1 — a’g’) +«°q = 0, (1) 


known as the “van der Pol differential equation”. This differential equation cannot be 
solved in closed form. Properties of its solutions have been studied extensively, however, 
by various means and from different points of view. Sparked by a paper by Levinson 


and Smith [1], dealing with 
q+99,0)7 +h = 9, (2) 


. theory for this rather general class of differential equations was developed [2]. Because 
the assumptions made in that theory in regard to the coefficients of the differential 
equation are rather weak (in order to keep the conclusions general), no very specific 
results can be expected from it. Hence it seems desirable to investigate special differ- 
ential equations or classes of differential equations, describing self sustained oscillations, 
which can be treated more fully; preferably differential equations for which closed form 
solutions, at least for the first integral, can be provided, if only in form of quadratures. 

One of the authors [3] has drawn attention to two such differential equations, which 
may be called “modified” or “‘associated’”’ van der Pol equations. They are 


q° — (sen q')(e/2)q°(1 — a’) + f(g) = 0 (3) 
ant 

— (sgn g°’)(x’s)(1 — 6’q°*) + «’f(q) = 0 (s > 0). (4) 
Equation (3) [describing a first special class of self sustained oscillations] has been 


treated by the authors in a recent paper [4], and a more general form of that equation, 


q° — (sgn q )(e 2)q°*h(q) + « f(q) = 0, (5) 


vhere h(q) denotes any suitable even function, has been subsequently considered [5]. 
Equation (4) [deseribing a second special class] is the object of the present study*. This 
equation has the convenient feature of allowing a discussion of its first integral by 
ther simple means. 

The second term in (4), which governs the energy input and output, is shown as a 

‘Received January 6, 1959; revised manuscript received May 19, 1959. The results presented in 
this paper were obtained in the course of a study supported by a grant from the National Science Foun- 
dation. 

Professor of Engineering Mechanics, Stanford University, Stanford, California. 

Professor of Mathematics, Ohio State University, Columbus, Ohio. 

‘Some of the results described here, mostly for the special case n = 1, can also be found in the 
recently published book by H. Kauderer [6]. 
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Fia. 1. Second term in Eq. (4) 
n = —(x*s)(sgn g')(1 — 6%q"*). 


function of 8q° in Fig. 1. In a way familiar from the discussion of the van der Pol Equation 
(1) it can be shown that the oscillations governed by Eq. (4) also tend towards a station- 
ary regime, a limit cycle. It is the purpose of this paper to study the maximum dis- 
placements (“amplitudes’’) Q of the oscillations governed by Eq. (4) in general and give 
particular attention to the amplitude Q* of the limit cycle. 

2. First integral; amplitude relationship. In Eq. (4), g denotes the dependent variable 
dots indicate derivatives with respect to time; x’, 8°, s are positive 


(displacement 
function f(g), which describes the restoring force, will be subjected to 


constants. The 
the limitations 
f(q 20 for g>0 


f(g) £0 for q<0 


f(q) #0 
so that only “genuine’’ restoring forces are admitted. Furthermore, we will start our 
discussions by requiring f to be an odd function, 

(-@ = {9 (7) 
This restriction, however, is imposed solely for the sake of simplicity of the presentation; 
it is not essential and will be dropped later (Sec. 4). 
3ecause of the factor (sgn q'), there are two different forms of the differential equa- 
tion, depending on the direction of the motion. However, if we reverse the direction of 
the coordinate g when the motion changes direction, by putting 


Gg=4q for gq <0 g 
q = —-q for q > 0, 
Eq. (4) assumes the single form 
q°+ «sil —@6(7)) + «f@ = 0, (9) 


with ¢ as the dependent variable. 
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Writing U = @°/x’ we find from (9) 

dU ee 

ar™ — B«U) + 2f(q) = 0, (10) 
a first order differential equation for U as a function of g. By putting 

aie U 222 , 
2 = G/s: y=; = B«s; “i = ¢(z), (11) 
S 8 

we may give differential Equation (10) a non-dimensional form, 


oY _ ry = —2[1 + o(@)]. (12) 


For the initial condition g = Q, > 0, when q’ = 0, equivalently z = X, for y = 0, 
the solution of (12) can be stated as 


y = 2e[1(X,) — I(a)] (13a) 
with I(x) denoting 


(a) = [ e™[1 + oO] ab. (13b) 


0 


Equation (13a) holds until the motion comes to its next stop (y = 0) atz = X, < 0. 
Hence 
I(X.) = 1(X,) (14) 
provides the relationship between two consecutive maximum displacements Q; > 0 
and Q, < 0 expressed by their non-dimensional representatives X, and X, . We shall 
refer to (14) as the “amplitude relationship” (taking the liberty of using the term 
‘“‘amplitude”’ as synonymous with “maximum displacement’’). 
3. Amplitude relationships for special types of restoring forces. We start our dis- 


cussion by assuming the function f(g) to be proportional to a pure (odd) power n = 
2m + 1 of q; 


f(q) = uw" ‘gq’,  (nodd). (15) 
The function ¢(zx), then, turns into 
d(x) = 7,0" with y, = (us)"™”; (16) 
the expression J(x) in (13b) becomes 


I(x) = 5 le~*— 1] + v | e ~ ede, (17) 


Carrying out the integration and transferring to the left hand side all terms independent 


of x, we find 


1} _ (2,)""" (2n)" _ -me| (Ne : e)") 
Yn | n! I(x) + n! al | sete 7y,n! a p> y! ° (18) 
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Using the symbols z = 2Ax and Z; = 2\X; we may abbreviate the right hand side of 
(18) by G,(z), where 





G,@) = e’F,(@), (19) 
FQ = a+ 5, (20a) 

and 
italia in (20b) 


Now, the amplitude relationship (14) may be written as 


G,(Z2) = G,(Z,). (21) 


Because F,(z), Eq. (20a), is a polynomial of odd degree having real coefficients, 
F,,(z) possesses at least one real zero, z) ; and, because the coefficients of that polynomial 
are positive, z) has to be negative. Furthermore we note that 


G,(z) > 0 for 220 

G,(z) — 0 as 2-70 (22) 

G,(@)—-7~—-37 as z7—o%, 
There exists only one abscissa z, for which G’(z) = 0; it is 

z, = —2nr/y,’". 

Since G’’(z,) < 0 the function G,(z) has a maximum at z, . From this and (22) it follows 
that 

2<s,, @ |e|>I|z, |. (23) 


Figure 2 shows a sketch of the general behavior of G,(z). 
A special consideration leads to the estimate for | 2 |, 


|Z | < max (n,a,) = M,. (24) 
The reasoning runs as follows. From Eq. (20a) we obtain for negative arguments z = —f 


(with ¢ > 0) 
eee es eee r( -§) Zs — ( ~2 i 
dil r+ h(i +e Ble ye | Balt a 


and we want to find a lower bound for the values ¢ that render F(¢) negative. If ¢ = n, 
the last term in (24a) vanishes, all other parentheses and, hence, {-terms are negative; 
therefore, for ¢ = n we find F(¢) < 0, provided a, < n. But if a, > n then F(¢) < 0 
for all ¢ = a,. 

Hence, we reach the following conclusions in regard to the amplitudes of the oscil- 
lation (illustrated by Fig. 2). 

(1) The amplitudes Q, , Q; --- (being proportional to X; and Z;,) are limited for all 


(finite odd) values of n in (15) and for any starting value Q, . In fact, for any 7 2 2 the 
inequality 
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Fig. 2. General behavior of function G(z). 


S 
ss. '3s/12\/8=3 (2 
Y: Q:| 5 5M. (25) 


holds, with M,, being given by (24). 
(2) For two oscillations starting from values Q, and Q, respectively, where Q, > Q, > 0 


we find 
1Q:1>1Q | forall 52 2. (26) 
(3) For any (finite) value of n there exists a uniquely determined limit amplitude Q*. 
The limit amplitude Q* obeys the inequality (25). More specifically, Q* = (s/2d)z* is deter- 
mined by the root z* of the transcendental equation 
G,(—z) = G,(2), (27) 


which, because of Eqs. (19) and (20), is equivalent to 


(n—1)/2 re (n-1)/2 2°? 
tanh z = | p> ap |/o a > |, (28) 


where for n = 1 the denominator 1s a, . 
Equation (28) represents a relationship between a, and z* as a function of a, (i.e. 
of the parameters of the system); we profitably use (28) for establishing the inverse 


function a,,(z*). From (28) and (20b) it follows that 2* is a monotone increasing function 


of X. 
Now, let us consider special values for n. The case n = 1, because of y, = 1, leads to 
G,(z) = e*[a, + 2] (29a) 
with 
a, = 1+ QZ. (29b) 
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Fic. 3. Functions G,(z) and G,( —z) leading to limit value 2*, 


Figure 3 applies to this case with X = 1/2. From the figure we find the sequences of 
amplitudes. If Z; > 2*, the sequence decreases; if Z, < 2*, the sequence increases to the 
limit evcle value 2*. From (28) we obtain 

tanh z ] 
— = (30) 
2 a, 
as the 2* — a, — relationship which determines the limit amplitude z*. The plot for 
Eq. (30) is shown in‘Fig. t.€For the example 2\ = 1 we find 
‘ >= (30: 
tanh z 2 30a) 





5 











Fia. 4. a; — 2* — relationship. 
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Since tanh z < 1 we have z* < 2, and the numerical value is z* = 1.915, equivalently 
zz” 1.915 or Q* = 1.915 s for the case considered. 
Figure 3 is plotted for the parameter value 2\ = 1. The amplitudes of the two sample 


sequences are 


Z;° =5, Z," = —1.99, Z;' = 1.9% 
Zz, = 05, Z, = —-1.73, Z;" = 1.90. 


For large values of the parameter 2, Eq. (30) may be replaced by 





2* wa; (31) 
from this there follows 
x*w14+2 or Qe a(1 + 1). (32) 
2X 2d 
Next, we consider the case n = 3. It leads to 
(2n)* 
a, = 1+ = 
and to 
G(z) = (a, +z+ Z + =). (33) 
From it the transcendental equation 
tanh z = 6z + 2° (34) 


6a, + 32° 





z Patt 


ye 









































[ 
0+ ’ : : - T T T T 
0 5 - |= 


Fia. 5. a3 — 2* — relationship. 
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for the limit value z* is derived. The plot for the z2* — a; — relationship is shown in Fig. 
5. For 2\ = 1, yz; = 1, we find a; = 7/6; hence Eq. (34) becomes 





> 3 
6z + z (35) 
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Fia. 6. Functions G;(z) and G;( —z) leading to limit value z*. 


= 


The plots of G;(z) and G; (— z) are shown in Fig. 6 for the parameter values, 2A = 1 


and y; = 1. From these plots e.g., the sequences 
zy = & , Z; = —1.%6, Z; = 1.55, 
Z> =0.5, Z? =-146, ZS = 1.52, 


for the dimensionless amplitudes can be obtained, both leading to the limit value 2*, 
given by the root of Eq. (35). That root is z* = 1.532, in agreement with the sequences 
shown above. 
Other special cases of odd values n in Eq. (15) can be treated in the same manner. 
4. Generalizations. Cases, where f(q) is not a single odd power of g, but a polynomial 
in odd powers, can be treated along the same lines. For instance, if 


fq = wg" + aa”, 


~ n,! z No! z 
mM ©14+4.— 257th =" (36) 
(2r)™ x v! (2r)"* x v! 
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Specifically, for n, = 1 and n, = 3, because of y, = 1, Eq. (36) turns into 
Pe) =1+204+29+4+ Sti (i+2+5+25) (37) 
oe 7 ae 2!" 3! 


which may be replaced by 
F(2) = (2d)*° + (2d)? + 6ys + 2((2d)? + Gys) + 3ys2” + 22”. (38) 


For the special values 2\ = 1, and y;, = 1, for instance, the equation (27) reads 
I ’ 1 


- 3 
iz +z 
tanh z = —— 3: 39 
8 + 3,7 ( ) 
This transcendental equation has the root z* = 1.282 leading to the limit amplitude 


()* 282 s. 

If we drop the restriction, imposed by Eq. (6b), on the function f(q) to be an odd 
function, and admit non-odd functions f = fy + f, [composed of odd, fp , and even, 
f, , parts] the function f(q) in (9) will differ for the swings in the two directions. We find 


(M@ =f(9 +f(q for gq <0 (40a) 

fq) = f(g) — fq) for g >O0. (40b) 

Hence, the swings with negative velocities, gq < 0, and those with positive velocities, 
q > 0, will have to be treated separately. The sequences | Z, |, | Z, | --- will lead to a 
value 2* developed from expression (40a) in the differential equation; the sequences 
Z,|,| Z, | --+ will lead to a value z2* developed from expression (40b) in the differential 


equation. The limit amplitudes Q* on the left hand side, derived from z¥ , and those on 
the right hand side, derived from z* , therefore, in general, will have different absolute 


values. 
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engineer, however, is that the author has included sufficient sensible examples and comments to give 
motivation to the mathematical exposition. 

Probably every reader of a book of this sort would make his own selections of material to be covered 
and the reviewer feels that a discussion of the normal distribution function which is essential to an 
understanding of statistical communication problems and of the unit impulse which is probably more 
used in this country than the unit step would be valuable to many readers. In general, this book will 
be most helpful to practicing engineers desiring to increase their mathematical competence and to 
advanced students needing a convenient reference covering most of the mathematics required in graduate 
engineering courses. The applied mathematician may also find the sections interpreting the mathe- 


matics in terms of engineering problems interesting. 
B. HazE.tTine 


Operational mathematics. By Ruel V. Churchill. Second edition. McGraw-Hill Book 
Co., Inc., New York, Toronto, London, 1958. ix + 337 pp. $7.00. 


Applied mathematicians will not need to be reminded that R. V. Churchill’s “Modern Operational 
Mathematics in Engineering’’, first published in 1944, was the first really accessible book on operational 
methods and the underlying theory of the Laplace transform to become available. The present book 
is a thorough revision of that earlier work. The chapter headings remain unaltered but much new material 
has been added especially to the last chapter. In the first two chapters, for instance, some of the oper- 
ational rules are derived by new methods, and in the third chapter problems involving servo-mechanisms 
and additional ones on circuit theory are introduced. The central six chapters of the book remain virtu- 
ally unaltered. It is in the last chapter that the major additions are made—the introduction of “gener- 
alized’’ Fourier transforms. Anyone who knows Professor Churchill’s writings will not need to be advised 
of the lucid manner in which he presents his material. The publishers have taken the opportunity to 
produce the new edition in a much more attractive way than the original and the result is a book which 
is a delight to use. It can be unreservedly recommended to anyone beginning the study of the theory 
of the Laplace transform and the operational calculus based on it and to anyone entrusted with its 


teaching. 
I. N. SNEDDON 


An introduction to plasticity. By William Prager. Addison-Wesley Publishing Co., Inc., 
Reading, Mass., 1959. viii + 148 pp. $9.50. 


The contents of this book are based on a series of lectures delivered by the author at the Federal 
Polytechnic Institute of Zurich in 1954. The Swiss edition of the book (in German) was published in 
1955 and the present volume under review represents more than a mere translation of the previous 
edition. Not only is this book broader in scope but it contains further exposition of the topics previously 
treated, a typical example being a more detailed treatment of the principle of virtual velocities. The 
references are brought up to date, an index added, and a large number of problems, totaling about 60, 
are included at the end of each chapter, thus making the book more suitable as a text. The material 
covered is very clearly presented, and the examples are carefully chosen. 

As in the previous edition, the book contains four chapters entitled ‘‘Mechanical Behavior of Plastic 
Solids’, ‘‘The Mechanical Behavior of Plastic Structures’. ‘Limit Analysis and Design’’, and “Finite 
Plastic Deformation’’. Chapters 3 and 4 contain several additional topics of recent origin not included 
in the previous edition, for example, a discussion of uniqueness, a discussion of two-parameter states 
of loading, as well as recent developments in connection with estimate of deflection at load-carrying 
capacity, and the theory of locking materials. 

By limiting the scope of the book to selected topics from the theory of perfectly plastic solids, the 
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LARGE AMPLITUDE OSCILLATIONS OF A TUBE OF INCOMPRESSIBLE 
ELASTIC MATERIAL* 


BY 
JAMES K. KNOWLES 
California Institute of Technology 


1. Introduction. In recent years a number of problems have been solved in the non- 
linear theory of elasticity for incompressible bodies (see [1] and [2] for discussion and 
references). The work of Rivlin and others has shown that the assumption of incom- 
pressibility makes possible the solution of these problems without any special assump- 
tions as to the form of the strain energy function characteristic of the material. Because 
of the severe non-linearity involved, the special problems considered have been static 
ones and the methods of attack have been essentially inverse. The question of waves in 
incompressible materials has been considered by Ericksen [7]. 

The present paper treats the dynamic problem of axially symmetric oscillations of 
an infinitely long circular cylindrical tube of incompressible elastic material. The sym- 
metry of the motion and the condition of incompressibility combine to permit the 
reduction of the problem to one to which the methods of the theory of non-linear vibra- 
tions of single-degree-of-freedom systems can be applied. The approach is semi-inverse 
in the sense that the determination of all physical quantities is made to depend on one 
unknown function which satisfies a non-linear ordinary differential equation. 

The specific problem considered here is that of arbitrary amplitude free oscillations 
of the tube when it is set in motion under given initial conditions and subject to no 
surface pressures. It is shown that periodic motions are possible for a large class of strain 
energy functions, and an expression for the period of oscillation is given in terms of the 
strain energy function. The results are shown to take a particularly simple form in the 
limiting case of a thin tube. Explicit results are given for the strain energy corresponding 
to a material of Mooney-Rivlin type. 

The static problem of the symmetrical deformation of a tube of incompressible 
material by uniform internal and external pressures is a special case of a more general 
problem solved by Rivlin [3] and discussed in [1]. It may be remarked that results similar 
to those contained in this paper can be obtained for the radial oscillations of a spherical 
shell of incompressible material. The corresponding static problem of the symmetrical 
expansion of a spherical shell under pressure was solved by Green and Shield [4] and is 
also discussed in [1]. 

2. Formulation of the problem. Consider an infinitely long circular cylindrical tube 
of homogeneous, isotropic, elastic, incompressible material which in its undeformed 
state has inner radius r, and outer radius r, . A point in the tube which at time ¢ has 
cylindrical coordinates R, 6, z is assumed to have been at the point r, 0, z in the un- 
deformed state. The motion is thus completely described by the function R = R(r, 2). 


*Received March 6, 1959. 
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If R,(t), R,(t) denote respectively the inner and outer radii at time ¢, the incom- 
pressibility condition (see [1]) may be written 


(2.1) 


ot 
4, 


lA 


R? —RI =r —-7r,, t= 0; , r fs. 
The motion is completely determined if R,(¢) is known. 

Since the formulation of the present problem differs from that of the corresponding 
static problem only by the inclusion of the inertia terms in the differential equations of 
motion, the statement of the complete system of equations of incompressible elasticity 
and its detailed reduction are omitted and the reader is referred to [1], pp. 88-92. The 


notation used in this section will be similar to that of [1] 


Let 
J=r/R 2.2) 
and temporarily regard Q and ¢ as independent variables. 
For an incompressible material the strain energy W per unit undeformed volume 


is a function of two of the three principal strain invariants J, , J2,/;: W = W(/,, J.). 
For the symmetrical deformation (RQ, 6, z) — (R, 6, z) considered here it is found that 
L=1=1+@+4+1/Q’. (2.3) 


Incompressibility requires that J; = 1, a relation which can be shown to be equivalent 


to (2.1). 
When the stress tensor is calculated in terms of W and the arbitrary hydrostatic 


pressure p, it is found that all shear stresses vanish identically. When the remaining 
stresses are introduced into the radial equation of motion, there follows the differential 


equation 


f) wy ow 2 (a + a 
— 2Q° — + 21+ Q’ _ + ¢ (2.: 
0Q E + 2Q ™ 2) oF | Q ) al ) 
__RQ_ aR 
a ] oo ot’ 
where p is the density. The boundary conditions at the inner and outer surfaces take 
for form 
Q=: pt2Q a + 21 + Q) > a = —P,(t), k=1,2. (2.5) 
LE Ol; 
Here P,(¢) and P,(¢) are the inner and outer viata pressures, respectively.* 
The incompressibility condition (2.1) is now used to compute the acceleration 
0°R/at’ in terms of the acceleration d’R,/dé’ of particles on the inner surface. Thus 


RQ aR Q aR, | 2 Q R? (2). dist 
a ho oe + =a (2.6) 
a” dt” i-¢ -@ Ry — ri j\ dt ; 





i-@ wr i-@¢ 


Equation (2.6) is now introduced into the right hand side of (2.4), the resulting equation 





*The boundary conditions at the ‘“ends’’ z = + © of the shell are taken to be of the form of vanish- 
ing radial and circumferential shearing stresses and vanishing displacement in the axial direction, so 
that axially symmetric motion which is independent of z is possible. This of course requires the presence 


of a suitable axial stress component on the ends of the cylinder. 
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is integrated with respect to Q and the boundary condition (2.5) at Q, = r,/R, is imposed. 
It is found in this way that the hydrostatic pressure p must satisfy 


_ Q 2 r 
2s on OW , aw 1+Q (2m am 
Ft Gg tte. = wa a 


#R, QdQ_ 
Oe Luwt+¢ = 


r (ae) fp (2 5- Om i.) de — P, (8) 


The remaining boundary condition then provides the differential equation for R,(¢) 
Thus when (2.5) is imposed at Q. = r/R, and the integrations in (2.7) are carried out, 
the following differential equation results. 


E | (daa Lae “a E inde) __ par (4) 
Of gt 108 R hed in ‘i Ri re — ri + Ri\ de 
Oe cca aiens fd 1 =z 7) _ 
re 2] 7 ao dQ = P,(t) — Pt). — (2.8) 


In (2.8) the incompressibility condition has been used to eliminate R, = (r.°—r,°+R,")'”. 
The inner radius R, and the velocity dR,/dt are supposed given at time ¢ = 0. 
3. The basic differential equation. The differential equation (2.8) is now rewritten 
in terms of the following quantities: 
z(t) = R,()/, | 
w= (i/r) - 1), 


bey Ah Ae (aw am) ] 
mee x, | Q al r O1,/; =I,=1+Q?+1/Q? seis ea 


In this notation (2.8) becomes a non-linear second order differential equation for the 














(3.1) 


dimensionless inner radius c(t): 


si (1 +4) £2 4 [ioe (1 44) - to (4 + te, = BO = PLO 
x log ( =) ap a | to 1+ 2 n° | di + f(x, nu) = rp/2 (3.3) 


In the sequel (3.3) will be considered only in the case of free oscillations [i.e., P,(é) = 
P,(t) = 0} and subject to the initial conditions 
x(0) = Xo ’ 


a (0) = (3.4) 


The function f(x, u) which is defined in (3.2) and which appears in the differential 
equation (3.3) can be rewritten as follows. Let 
2 


W(u) = a WU, , 1e)s,0t001+0st/s « (3.5) 


ip 








74 JAMES Kk. KNOWLES [Vol. XVIII, No. 1 


Then 


Ws 2 (aah, Wat) 2(, 1), aM) gg 
du -—s rip \ Al, du ol, du/ rip u/\ al, PS sie 


If this fact is used in conjunction with the change of variables Q = u~’” in the integral 


’ “a d - 
f(r, uw) = / (u — 1)7' av du. (3.7) 
(ptz?)/ (etl) du 


4. Free oscillations. With the notation v = dx/dt, d’x/dt? = v dv/dz, it is possible 
to write the differential equation (3.3) in the form 


(3.2), there follows 


free ee 
a [sux log (1 + uw/2’)] + af(z, uw) = O (4.1) 
ax 
and thus to obtain the first integral 
Ay?” log (1 + w/a’) + F(a, vw) = C. (4.2) 
The constant C = 4,25 log (1 + w/2>) + F(xo , uw) depends on the initial values x, and 
v) , and the function F(z, u) is given by 
~ ite . r nl dW, 
F(z, wu) = | Ef(é, wu) dé = | é | (u — 1) a du dé. (4.3) 
1 J1 J (w+? p+l) au 


If the order of integration in the repeated integral in (4.3) is reversed, it is possible to 
carry out one of the integrations and obtain 
F(x, w) = 3(x° — 1) | (u — 1) °W,(u) du, (4.4) 
J (u+z?)/ (ul 
a form which will be convenient later. In obtaining (4.4) it has been assumed that the 
strain energy vanishes in the undeformed state, so that 


W(1) = 0. (4.5) 
It is well known from the theory of vibrations (see, for example, [5]) that the motion 
x(t) is periodic if and only if the “energy curves’ (4.2) are closed curves in the x — v 


plane with a finite period £dx/v. We proceed to consider the circumstances under which 
periodic motions are possible. 
It is now assumed* that 


ow ow) 
— + - > 0 (4.6) 
(a : ol, 1,=I, ) : 
so that according to (3.6) dW,/du has the sign of u — 1 for u > O. From this and (3.7) 
and the fact that u > 0, it follows that 


a0 +f @<¢ <1 
=0 if rz=1 ‘ (4.7) 


f(a, pw) | 
i>0 #4 2z¢>1 





*This restriction on W is implied by inequalities proposed by Truesdell [2], p. 182, and is identical 
with an inequality obtained by Rivlin and Saunders [6]. 
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Consequently the function F(x, uw) defined in (4.3) decreases monotonically in x for 
0 < x < 1, vanishes at x = 1 and increases monotonically for x > 1. 

For real motions x(¢) it is always true that F(x, uw) < C, so that for every x such that 
F(x, uw) < C, Eq. (4.2) determines two values of the velocity v: 


1/2 


2C — 2F(zx, yu) | 

—s 

| 

x log (1 + 4,)| 

\ a} 

lo show that the energy curves (4.2) are closed, it is only necessary to show that 
there are exactly two values of x such that v = 0 in (4.2); i.e., such that 


F(a, uw) = C. (4.9) 


(4.8) 





v= +t 


decreases monotonically in x for 0 < x < 1, there will be exactly one 
positive root « < 1 of (4.9) for every positive C unless lim x — 0 + F(a, yu) is finite. If 
this limit is finite, values of C which exceed it will fail to determine a root x < 1 of (4.9). 
Similarly since F(a, u) is monotone increasing for x > 1, there is precisely one root 

> 1 of (4.9) unless lim x — © F(z, yz) is finite. If either of these two limits is finite 
periodic motions can exist only for sufficiently small C. 

The conditions that guarantee that the energy curves be closed (and in fact that the 
are thus that the integral (4.3) defining F(x, u) be unbounded as 
- and as x — o. To translate these conditions into requirements on the strain 
energy W, , it is convenient to consider F given in the form (4.4). From (4.4) it can be 
shown directly that lim « — © F(x, yw) will be infinite if only W (uw) is unbounded as 
lor the case lim « ~ 0 + F(z, ») it is observed that by the definition (3.5), 


Since F(a, p 


motio! be periodic 


> 


7 


W(1/u) = Welw). (4.10) 


[f the change of variables u — 1/u is then made in (4.4), it is found that 


1—27 s'” :' Poe 
Fi(z,nu) =- ——— [ (u — 1) “W,(w du. 
* u 1)/ 


From this it follows that lim 2 ~ 0 + F(z, nu) = @© if the integral 


a x 


| (u — 1)-°W,(u) du 
J1+1/p 
diverges. This will be the case if W,(u) increases at least as fast as u as u — ©. Hence 
the energy curves in general will be closed if the strain energy is such that W,(u) becomes 
infinite at least as fast as u as u— ©, 
The two roots x = a < landx = b > 1 of Eq. (4.9) are the minimum and maximum 
amplitudes, respectively, of the oscillating inner surface of the cylinder. 


The period of the motion is given by 


- : > dr ae d fs log (1 + w/a’) |'” 
T=2] “= 2/ Rese dx. (4.11) 


Ja “a 


Since C — F(x, ») can vanish only to the first order at x = a and x = b, and since 2° log 
(1 + w/z’) is bounded, T is finite. 
5. The thin shell. A limiting case which is of interest is that for which up — 0. 
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Since u = (r2/r;)” — 1, the case of small » corresponds to the thin shell. The question 
specifically considered here is that of the limiting form of the maximum and minimum 


amplitudes of oscillation and of the period as np — 0. 
Equation (4.9), which determines the amplitudes x = a (u) < 1 andz = b(u) > 1 
may be rewritten as follows when C is expressed in terms of the initial values and after 


division of both sides by u. 
F(x. u) 1 1? log a + p/2a) 4 F(xo _, uw). (5.1) 


) 


_ a ae ee 
rm 2 m rm 


From the representation (4.4) for F(z, «) it is easily shown that for every finite positive x 





. Fla, p) ere ae 
lim ~-= > W (2°). (5.2) 
REE 7 2 
Also 
. 1 .. log (1 + w/z) l > 7 

lim ——— @ Bh (5.3) 

oi ae 7 y 
a and 


If a, and b, denote the limiting values as » — 0 of the amplitudes a(yz) and b(u), 
bp satisfy the limiting form of (5.1) as u — 0. 

W(x’) = vo + W2(2%). (5.4) 
From the property W,(1/u) = W,(u) it follows that the maximum and minimum limiting 
amplitudes a, and by are related by 


b, = 3 . (5.5) 
Ao 


0 


Now consider the expression (4.11) for period 7. In the limit » — 0 the results (5.2), 
(5.3) and (5.5) show that 


1/ao 
T= limT=2[ — i + Wola) — Wola)? de 


is the limiting value of the period of the motion as u — 0. 
It may be remarked that while the present discussion has obtained only the limiting 


values of the amplitudes and the period as u — 0, it is possible by further calculation to 
obtain the correction terms of order yu, say, in the formal expansions of these quantities 


for small y. 
6. The Mooney-Rivlin material. It is interesting to examine the results of Sec.4 


and 5 for the case of the Mooney-Rivlin material (see [1], p. 76) for which the strain 


energy function is given by 
W(, , J2) = 4a(J, — 3) + 3617, — 3), (6.1) 


where a and £ are positive constants. 
According to (3.5) the function W,(u) assumes the form 


W,(u) = K(u ++ ~ 2), (6.2) 
i 


where 


K = (a + £8) 1; p. 


When the integrations in (4.4) are carried out it is found that 
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. oe 1 + p/2” 
F(x, w) = 34K(1 — 2’) lo (Lt), 6.3 


For this special case F(x, u) is unbounded as x — 0 and as z — ©, so that periodic 


motions exist for all values of the initial conditions. 
In the limiting case of the thin shell the condition (5.4) for the amplitudes becomes 


2g] _%, 21 
et+a=Rtnmt a (6.4) 


The root a, < 1 of this equation is readily found to be 


= E 4% 1 a 1 (2, + 4 4 1 4 2x0 Quo 2) "]" (6.5) 
M~loK 2 Taf 2\ki'** AT K * Kz a 








while the maximum amplitude is given by bb = 1/ao. 
For the initial conditions 7, = 1, vo # 0, corresponding to the situation in which 


the tube is set in motion by imparting to it an initial velocity in its undeformed state, 


Eq. (6.5) becomes 





r 1 vs 4v?) “yl 
aA = E + -1(3+#) 
2K 2 K (6.6) 
b, = I 
ao 


On the other hand for the initial conditions x) ¥ 1, 1). = 0, so that the tube is released 
from rest in an expanded or contracted state, Eq. (6.5) shows that the inner radius of 
the tube oscillates between its initial value and the reciprocal of that value. 

The expression (5.6) for the limiting value of the period takes the form 


- ee Vee fy? , 1 : 1\7'”? 
T, = OE [ —+uytao ore dz, 
xo x 


Jac K 
where ap is given by (6.5). When this integral is evaluated it is found that 
2 1/2 
T» = a = ( 47 ) . (6.7 
"\a + 8 


Thus in the limit « — 0 the period of oscillation does not depend on the amplitude or 
any other characteristics of the motion for the Mooney-Rivlin material. 
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author has been able to lead the reader to the forefront of the current state of research on these topics 
without requiring prior familiarity with the subject. This has been achieved by presenting the basic 
theory in a condensed manner in the first 18 pages of Chapter 1 through efficient use of mechanical 
models. Those desiring a detailed account of the development of the subject based on the many experi- 
mental investigations and theoretical contributions, beginning with the work of Tresca in 1868, should 
refer themselves to the several other available books in plasticity. 

The chapter on limit analysis and design, comprising half of the book, (63 pages) is the most valuable 
part. Fifty references (mostly from works of the last decade) are cited, covering the fundamental theorems 
with application to frames, circular plates, and cylindrical shells. 

The reviewer is confident that the author has achieved his hope, stated in the Preface, that the 
“book will usefully supplement the traditional texts on strength of materials and theory of structures 
which are primarily concerned with the elastic behavior” 

P. M. NaGuoi 


Principles of optics. By Max Born and Emil Wolf. Pergamon Press, New York, London, 
Paris, Los Angeles, 1959. xxvi + 808 pp. $17.50. 


“Principles of Optics’’ is a new book; it is not simply a translation and revision of the famous “Optik’’ 
of 1933. This treatise is a model of completeness and thoroughness in classical optics, both geometrical 
and physical. It is a tremendous piece of work, and so well done that it is certain to be an essential 
part of the library of anyone concerned with optical phenomena. 

The book contains more than 800 pages with 14 chapters and 9 appendices. The first two chapters 
are an excellent review of electromagnetic theory with special attention to the reflection refraction 
and transmission of plane waves. Elementary dispersion theory is also dealt with at this stage. The third 
chapter, Foundations of Geometrical Optics, shows how geometrical optics is obtained from Maxwell’s 
equations as a limiting case for short wavelengths—the eikonal equation. The mathematics of geometri- 
cal or ray optics is introduced. Chapter 4 deals in detail with the geometrical theory of optical imaging 
using the various characteristic functions of Hamilton’s methods. Chapter 5 is concerned with the 
geometrical theory of aberrations. Chapter 6 is entitled Image Forming Instruments and deals with 
the human eye, the camera, the refracting telescope, the reflecting telescope, and the microscope. Chapter 
7 is a chapter of more than 100 pages on Elements of the Theory of Interference and Interferometers. 
Chapters 8, 9, 10, and 12 are’ devoted mainly to diffraction theory. Chapter 8, Elements of Diffraction 
Theory, is the best discussion of classical diffraction theory that the reviewer has seen in any book. 
Chapter 9, The Diffraction Theory of Aberrations, Chapter 10, Interference and Diffraction with Par- 
tially Coherent Light—this is one of the more sophisticated sections of Optics, and it is dealt with 
rather thoroughly here. The theory of partial coherence deals with correlation functions and time 
averaged intensities—the quantities one measures in experiment. This section of the book should be of 
interest to mathematicians as well as physicists. The same can be said of Chapter 11 entitled Regions 
Diffraction Theory. This chapter is confined mainly to the half-plane problem of Sommerfeld. It is 
here that the reviewer feels that the authors should have discussed the newer integral equation and 
variational techniques that lead to solutions of diffraction problems in closed form. This could have 
been done spacewise at the expense of several of the appendices, especially Appendix I on the calculus 
of variations, as well as Chapter 12 entitled Diffraction of Light by Ultrasonic Waves. 

The remaining part of the book is devoted to excellent discussions of Optics of Metals (Chapter 13) 
and Optics of Crystals (Chapter 14), The discussion of the Optics of Metals in Chapter 13 is the classical 
version and one might have asked for a brief discussion of the importance of quantum mechanics in 
this area. Altogether the reviewer feels that this is by far the best advanced text in Optics now available. 


Roun TRUELL 
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THE CALCULATION OF HEAT FLOW IN MELTING SOLIDS 


BY 
MARK LOTKIN 
AVCO Research and Development, Wilmington, Mass. 


I. Introduction. The numerical integration of certain heat conduction equations 
was treated by the author in [1]. The case considered there dealt with the transfer of 
heat in a one-dimensional finite slab composed of materials possessing variable thermal 
properties, subject to appropriate initial and boundary conditions. 

The purpose of this note is to extend the previous discussion to include the phe- 
nomenon of melting. In carrying out this extension it turns out to be advisable to modify 
the numerical technique of integration used in [1] by going to unequal subdivisions 
in both space and time variables. Near the melting surface, where the changes of tem- 
perature may be very pronounced, it is advantageous to choose rather small integration 
steps. Farther away from the melting surface, on the other hand, where the temperature 
may change rather slowly, a more economic treatment of the integration process may 
be achieved by taking relatively large intervals in the independent variables. 

In the following sections, then, the physical situation is first defined, and the govern- 
ing mathematical equations are stated. These differential equations are next replaced 
by analogous difference equations, and the resulting truncation terms are noted. It turns 
out that the stability characteristics of the method are quite similar to those established 
in [1]. 

As numerous calculations carried out on the IBM 704 machine by means of a program 
based on some of the equations discussed here have shown, the procedure derived in 
this note leads to quite satisfactory results. 

II. System of differential equations. Let us consider, then, the one-dimensional 
conduction of heat in a solid slab consisting of one material, and extending from x = 0 
to x = a, with the surface at x = a to be insulated. Through the face initially located 
at « = 0 let heat enter at the rate g(t). Let it further be assumed that the amount of 
heat flowing through the frontal face is at times sufficiently large to raise the frontal 
temperature to the melting temperature u,, of the material, thus actually setting off 
the melting process, and continuing it until the temperature drops temporarily below 
Um . If, upon melting, the liquid material is assumed to be removed instantaneously 
the “frontal’’ face of the solid is seen to move forward, to a position z = s(t) at time t. 
Of interest then is the temperature u = u(z, t), and the amount s(¢) melted. 

The basic equations governing the melting problem described above are: 


Ou te] Ou 
~-5(,5 <2< 0 1 
Pat ax (i ) Wwssse, +26, (1) 
where c, p, k denote known functions of the temperature uw. 
Initially, 


u(x, 0) = f(x) <u, 0<2z <a, i = 0. (2) 
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At the frontal face 


—k(du/dx) = g(t) — pFs'|* (3a) 
with s(t) = 0 for u(s, t) <u,¢r = s(2), t>0. (3b) 
s(t) > O for u(s, t) = u,, (3c) 

At the back face, 
(du/dx) = 0, x=a, t> 0. (4) 


As pointed out by Landau [2], the numerical treatment of the problem is simplified 
somewhat by the removal of the moving boundary, by means of the transformation 


& = (x — s)/(a — 8). (5) 
Clearly above equations then become. 


u(z, t) = U€, 0), 











UE, 0) = ft) <u, O<E<1 t= 0, (7) 
with 
“> or (22) ao = oe (8a) 
s(t)=0 for U0,) <ulg~=0, t>0 (Sb) 
s(t) >0 for U(0, t) = unJ (8c) 
(a@U/at) =0 €=1, €t>O. (9) 


III. Definitions and assumptions. Let us assume that n + 1 values &; ,7 = 0,1,---+n 
are given in0 < é < 1, with & = 0, &, = 1. Further, let Vié = Vi = & — 6-1, Vit = 
t; — t;-, fori = 1, 2, --- n,j = 1, 2, «++ 5 t = 0, and U(é, ,t;) = U;; , c(U,;) 
CG; , a(t;) = gq; , @ = k(ep)™", and r,, = a; Visit [((@ — 8;) Viet)’, with Vint = 
V,éfort =n. 


In the subsequent numerical treatment extensive use is made of the relationships: 


1 /d"w ‘— tm({D° "wh - o) sii 


™m: i=0 


which expresses the presumably existing derivatives of an arbitrary differentiable 
function w(z) evaluated at z = 2 in terms of the divided differences based on the points 
(z,, w,;),4 = 0,1,2,---: 


(Dw/Dz):o = (wi — Wo)/(@: — Zo) 
(D?w/ Dz’) oo = [(Dw/Dz)2, — (Dw/Dz),0]/(z2 — 20), ete.; (11) 


in Eq. (10) of” denotes the ith symmetric function of 79 , 71 , -** 5 Ti+m-1 5 T? = 2 — 2p, 
7 = max | z — z, |, and (d**’w/dz""') is to be evaluated at an intermediate value ¢. 


*The symbol s° denotes ds/dt. 
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Thus, 


aU 7" Uisr.ie1 — Ui js oa (Pe) 
(z ) - V iat Views De i+2,4 








(12) 
D’U , 5 
+ ViaVior + Vive) Hes + C, | Vel’, 
Dé i+3,¢ 
aU pe U5 541 =- Uy, a . 
(2 ) = Vial 7 Vin 2 i4+2,3 + C(VO) : (13) 
1(#U) —_ (DW) 
4(% ). _ (FE vag FT C2! VEL: (14) 
with 
|Vé| = max | V,é], Vi= max Vil. 
p=0,1,***,n 770,1,2,°°°,3 
Alternately, clearly also 
-(#@) - cts v (FE) +cuve, as) 
1 (#U - pr) ee 
5 (SF) 7 (7 iin (16) 


As to the nature of the function z = s(¢), it must clearly be continuous, and possess 
a non-negative derivative &(¢) except for a finite number of values t,; at which either 
$(¢;.) or &(t;-) vanishes. 

IV. Equation for the melting rate. For the following discussion it is convenient 
again to re-write Eq. [6] in the form: 


aU _ (a —s)° ou _1=t, 2) _k (uy . 
a” a | a—s° a k Vag]? (17) 
with k’ = dk/dvU. 


In deriving suitable difference equations let us consider first the equation (8) governing 
the heated surface ¢ = 0. At time t;,, , by Eq. (8a), 


Ress aU “ 
oe (U) = O54: — Po Fs; ., 
0,j3+1 


a — 8j+1 \ 08 





For the following considerations it is necessary to distinguish two cases, depending 
upon whether the surface temperature at ¢;,, is below or at the melting point u,, . 
i) Ug as < Ua 
Then &;.; = 0, s;4. = 8; . Eq. (8a) becomes 
—(dU/0)o, 541 = (a aes 8 )koidies a (a cont s)(k~’).@ Vi. 
Now replacing (U;)o,;+1 by Eq. (12), and (D?U/Dé’)... by Eq. (14), (Uge)oui+. by 
Eq. (17), (Uso,;41 by Eq. (18), there is obtained 


(1 + 270;)Uo. 541 — 2roiU 1,541 = Uo; + 2ro(a — 8;) Vikoidiar +R, (18) 
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with 
te Wl_U} == g Ugvt |] _ «WV. Vt ae 
m= aVet k FE — s)° k(a — 8) _| 3(a — s)° : Veee + a U(Vty. 


a) Uscsus = tes 
a (U)o,441 = 0. Eq. (8a) may be written as 


i kU. : 
F poj8j41 = Qi41 2 cane 5 (Us)o,:41 + vi 2 meer U, — Fs » | 


Again replacing (U;).,;,: as indicated above under (2), there results, after some calcu- 


lations, 


8341(U go, 541 + R; , (19) 





Fp, 8341 = q | > nage (Um sions U, i+1) + 


where 


‘Or { S ’ C 
R; = vou - + > = tk VV; U; «| + v4 a 4 - U; zai Fe'a, | - V1 es = U ft - 


Xa — s) a-—s a-s 


Since s;., = 8; + Vysits's41 + 0 ((V2)*), Eq. (19) may be replaced by 


; | Ko; , iV ; 
a * So E asl ag (Uy — 1,441) rh =~ g; (UO so, ul +- Ry’, (20) 


Fp Vi(a — 8;) 2a; 





with 


ie «m4 2S [o, re re | 
(a — s) ? 2 


In Eqs. (19) and (20) (Ué)o,;.; is to be evaluated by means of Eq. (12). Since the accurate 
determination of s° is essential for the closeness of approximation of the solution of 
the difference equations to the solution of the differential equations, the evaluation 
of (Ué),.;., by at least a four point formula seems advisable. 


Equation (19) may be expressed in the form 


Sia. = Ajar — Byala — 8;4:)°° + RP, (21) 
where 
TAO esos 1 
A j+1 = tres] Fou — 2 sis a Vid ics | > 0, 
“Qo; 
ee | ee | .  k vie |” 
esse ee RE eS 7 ‘sini’ ‘guietaammaaaoemae ( 
Bis V; Fp 2a > 0, 
b r a | 
R, = Rt ‘ - EvU | 
V4 6 4 


It is assumed here that (Ué)>.;.; < 0. 
Similarly, from Eq. (20) 
Sj41 = Aja: — Byii(a — 8;) + RE. (22) 
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The numerical solution of Eq. (21) may be accomplished by the iterative techniquef 


>” = A — Baa — o™)" 


of) = go, + (Visit/2(2; + =), (23) 
for n 0, 1,2, +--+, witho® =; and >> = de/dt. It is seen that 
n41 n) B V jail ( (n) lies 
Cg = oF = ee as. tan-th o ~ . 
Aa — o”™)(a — a”) “ 
Let us assume now that a — o” > b holds for all n; this is obviously a practical 
assumption. Then 
isae ; BViait n) — 
oe? — 6 | $e |e — oI, 
whence 
> ae! ae aw (0) | BV i4:tY : 
o = ¢ Ss (2 Yint) o —@ = (BY! | V j+ilZo p 
The sequence o’ is thus convergent provided | BV,;,,t/2b° | < 1, a condition which 
is found to be satisfied in practice. Thus, for example, for the values a = 1 inch, b = 
1/3) inch, V, = 1/100, u, = 5000, p = 150, ¢ = 1/2, k = 1/3000, F = 3000, V,.:¢ = 
1/4, and U, ;;; 1000 there results BVt/2b° ~ 1/100, so that the iteration will con- 


verge rapidly. If we denote lim,.... ¢;"" by o,; , then it is seen that 
¢,; — 8; | < Const V;t, for O< 7 < ju - 
V. System of difference equations. The development of the difference equation 
analogous to the differential equation (17) for interior points is somewhat more laborious. 
Applying Eq. (17) at the point &; , t;., , replacing 0°U/aé* by Eq. (16), 9U/dt by Eq. 


(13), and 8U/d& by the central divided difference approximation obtained from Eqs. 
(12) and (15), 


au Lk ay ee ie ‘Pe Ai 
— Af Ussrier ~ © iiot UC iijen ~ Ui-n.isi ’ e | 
(22) “8 a] V cas ¥ V Jee vel 


there is obtained, after some calculation, with pj = Visi/Vi , 
yp 1-ég .. | , 
— | 2r,,; Pe — 8 Viet [Winns 
| Titp, W2Za-—-s)V;i' 7 “are 


s 1 — é; P : 
we [ + 2p7,, — (i - Dp: ) - 3) V, 8; Vint fe are 





(24) 


al iy | + p, Ha — s,) V,P'* 7+1 i+1,j+1 


U.; + rs (’) Wiess — 0 — 00a — BU ral 


+ V ACs | Ve + C.(V.)]. 


For convenience in writing the subscript j + 1 has been deleted in some of the following equations. 








84 MARK LOTKIN [Vol. XVIII, No. 1 


Let us, finally, turn to the condition (9) expressing the insulation of the rear boundary. 
The difference equation representing this condition may be’ obtained either gpd 


(13), and (17), or by applying Eq. (24) to the case i = n. 


by means of Eqs. (12), (15), 
Ui-1,4, Vaei = Vo, there is vad ained 


In this case, since &, = 1, and assuming U,.,,; = 
immediately 


2rnjUnar,ja1 + (Ll + Qaj)U nie. = Uni + VelC, | VE| + CVO). (25) 


2r nj 


It follows, eile that the system of differential equations (1), (2), (3) may be approxi- 
mated by the following system of difference equations: 


(1 oo 2r Ze Ee, as 2r he == T 0; oh 2ro;(a —_ 8;) VY Bertie if T'o.¢41 <= Un (26) 








Dis = Ager — Byala — 0541)", (27) 
Oss = 0; + (Vaart/2)[2; + 2j%], n=0,1,2,---, (28a) 
eee = Um ’ (28b) 
= 2.57 21,441 + (1 + 2rni) Tn, 541 se ie (29) 
Di 1l—é ; “ 
ar aes ey 25 Viet Pe-rs01 + 
: + 2, 2a—a;))V; ° ” 
bli +o tt erie oe y. V...t I7 
I I a! ;;D; Pi Xa —a,;) V; 7+1 t,7+1 
(30) 
Di i-é : 
— | 2r,, —& i pty V afr il 
reste cena SO it) wl lade 
ge Whe’ ae = . 
‘ty = ( ‘) [Ti41.5 — 1 — pdT i; — piTi-1,;]. 
t k tj : 
Equation (30) applies toz = 1, 2, ---,n — 1. 
Equations (26), (29) and (30) revert to the corresponding Eqs. (10), (12), and (8) 
of [1] in case of non-melting; it is necessary only to replace aV,; by Az, Vt by At, a; 


For equal values of Vx and V¢ the equations become essentially identical with those 


of [3], with the exception of the relationship (19) and (20) for the rate of melting. 


VI. Integration procedure. The onset and stopping of melting presents some minor 


problems in the solution of the system (26) through (30) of difference equations. A possible 


approach is the following: 
(a) Initially, at 7 = 0, o = 0, 2 = 0. Eqs. (26), (29), and (30) represent n + 1 


linear equations in the n + 1 unknowns 7; , for 7 = 0, 1, --- , n. The matrix of the 
linear system is tri-diagonal, and may thus be treated by the method of solution of [1]. 
This step is repeated, for an arbitrary choice of time increments V/,,,¢, so long as 


Feds < By. 
(b) Eventually it may happen that 7; = T, (¢;) < u,, but T4541 = Tot; + Vi4it) > 
Then the solution of Eqs. (26), (29), and (30) is carried out for a progressively 


Un 
reduced sequence V;.,¢, until 75,44: = To(t;+:) satisfies | 7’ (t;4:) -un| < ¢6€ > 0 
prescribed. The values of T. +1 10r 4 be , n are now assumed to be correct. 


.. $— 
Placing in Eq. (27) n = 0,¢;3; = 0; p Th: = : » Visit = tar —t;, and (T'g)o,j41 
1 


evaluated from 79,541 = Um» 11,441 » T'2,4+1, 7'3.4+1 , by Eq. (12), there is obtained the 
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value of >-,") , and then o,°}} by Eq. (28a). This operation is repeated until convergence 
to limiting values o;,,; > 0 has been achieved. 

(c) The integration is now resumed, for the equations (28b), (29), and (30), where 
c; , >.; are replaced, respectively, by the values of o;,; and > ;,, obtained by (28a), 
(27). While V/;.,¢ is arbitrary, it should be restricted by the requirement that V;.:s < 
V2 or >>; Viait < (a — ¢;) V1. Then Eq. (27) is utilized again. 

This step is repeated so long as >. ;4, > 0. 

(d) Eventually it may happen that o;,, < 0. Equations (28b), (29), (30) are now 
solved with a sufficiently reduced value of Y7;.:¢ so as to produce >>;,, = 0, to within 
a prescribed tolerance. The next step then proceeds as in (a). 

VII. Applications. In testing the method described above one may conveniently 
use the case of the semi-infinite slab of a material which possesses constant thermal 
properties k, c, p. If it is further assumed that the heat input q, and the initial temper- 
ature distribution f(z) = wp are constant, then the solution, before melting starts, of 
the problem is given by [4] 


9 - ' és 
™ - rAd | . ; . ( tH ) bi ie | 72 | 
ulx, t) — = i (at BY — 1 . 
' to” be | (a xP \ dat a. SS 


For the time ¢,, at which melting starts, since ~(0, ¢,,) = u,, , there is obtained 


tn = mkep|(um — Uo)/2q]’. (31) 





Further, as Landau [2] pointed out, there exists a steady state solution of the system 
of differential equations (1), (2), (3) for the case of the semi-infinite slab; as t > @ 


ss. = g{[oF + cum — U)]}-, (32) 

s(t) os. = sat — kqo'(tm — Uo) (33) 
and, in general, 

© <4. S. < s(t) < s(t — t,,). (34) 


The results obtained by means of a program for the IBM 704 machine based on the 
system of difference equations discussed here have been quite satisfactory. In particular, 
Eq. (31), indicating the onset of melting, and Eq. (32) for the rate of melting, have been 


found to hold to high degrees of accuracy. 
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—NOTES— 


THE OSCILLATIONS OF A VISCOUS LIQUID DROP* 


By W. H. REID (Brown University) 


1. Introduction. In his discussion of the effect of viscosity on the small oscillations 
of a liquid globe, Lamb [1, pp. 639-641] observed that the results obtained in the limiting 
case of ‘small viscosity’ are independent of the nature of the forces which produce the 
tendency to the spherical form. When these forces are due to self-gravitation, the problem 
has been solved completely by Lamb [2] and Chandrasekhar [3] for arbitrary values of 
the viscosity, and the latter author has given detailed numerical results for the aperiodic 
modes of decay for / = 1, 2, 3 and 4, where / is the order of the spherical harmonic 
deformation considered. 

The question then arises whether or not these more general results for arbitrary 
values of the viscosity are also independent of the nature of the forces which produce the 
tendency to the spherical form. While no attempt will be made to answer this question 
in complete generality, it will be shown that when these forces are due to surface tension 
the results obtained are identical with those obtained by Lamb and Chandrasekhar for 
a self-gravitating globe. 

2. Statement of the problem. We consider a liquid drop which, in the undisturbed 
state, will be spherical under the action of surface tension forces. If the external pressure 
is taken to be zero, then the internal pressure will have the constant value 


p = 27;,/R, (1) 


where 7’, is the surface tension per unit length and RF is the radius of the sphere. To 
study the oscillations of this configuration we consider a deformation of the liquid 
surface of the form 


ile Ril + eY"(6, ¢)], (2) 


where Y7(@, ¢) = P7(@)e"'”* is a surface harmonic of the first kind and e < 1. In the 
absence of viscosity, assuming a time dependence « = «e*”’, the frequencies of oscillation 
are [1, p. 475] 


2 x 
gio = Ul — 1)(1+ 2), (3) 
pk 
and are the same for all values of m. 
When viscosity is included, the departures from the equilibrium state are then 


governed by the linearized equations of motion 


ou 6p 2 . / 
«fie —grad Ps yeur’u and divu=0. (4) 
c p 


*Received January 29, 1959. This work was sponsored by the Office of Naval Research under Con- 


tract Nonr 562(07). 
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On taking the divergence of this equation we have V*(ép/p) = 0; accordingly, the first 
order change in the pressure in the sphere is 


bp/p = €Pox' Y7, (5) 


where x = r/R and P, is a constant of integration that must be determined later from 
the boundary conditions. Equation (4) must now be solved for u as an inhomogeneous 
equation with 6p/p given by Eq. (5). In treating the viscous case, it is convenient to 
take the time dependence in the form e = ee “ so that the real part of o will be always 
positive. 

Since grad(ép/p) is a poloidal divergence-free vector, it can be represented in terms 
of a single defining scalar function II(z); in this representation its r-component is given by 


(crac 2p) = eo R ie . 6) 
p r d 


where 


I(r) = Mor’*’ and TM, = (l/o°R’)P, . (7) 


Since grad(ép/p) is a purely poloidal vector, it follows that u must also be purely poloidal 
and, accordingly, it too can be represented in terms of a single defining scalar function 
U(x); its r-eomponent can then be written in the form 
i= 
u, = eoR — i io (8) 
When defined in this manner, both II(x) and U(x) are dimensionless. The dynamical 
equation for U(x) which then follows from Eq. (4) is 


a... oe ee 
|S - +¢|o@ = ¢ne, (9) 


x 


where 


q= oR? /v. (10) 


3. Boundary conditions. At the deformed surface of the drop the radial component 
of the velocity must be equal to the velocity of the surface itself and this requirement 
leads to the kinematical boundary condition 


U=-1 at zr=1. (11) 
The other boundary conditions are obtained by considering the behavior of the viscous 


stresses at the surface. Thus, the requirement that the tangential stresses vanish at the 
surface leads to the single boundary condition 


d° 2 . 
E = 2 WED), =() at z= il. (12) 


da” x dx x 


These two boundary conditions are clearly independent of the nature of the forces 
which tend to maintain the spherical form, but the third boundary condition, which 
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follows from considering the normal component of the viscous stress, is not. This stress 
component is given by 


Ou 
: (1é 
or 13) 





—Prr = pt op — 2 
and the condition to be imposed at the surface is 
it 2 
—Drr asin n(2 +T x) ’ (14) 


where R, and R, are the principal radii of curvature, reckoned positive when directed 
inwards. For a surface deformed in the manner (2), we have [1, p. 475] 


1 1 1 . 
aqua a= = oun Poa 9 m / - 
R, + RR [2 + (2 — 1)(U + 2)eY¥%]. (15) 
The third boundary condition can therefore be written in the form 
i , d (4) . 
_ —— == — 2vo—\3 ‘ gm i. 
(1 1)(l + 2) oR Pox 2vo 7 \ at 2x 1 (16) 


The solution of Eq. (9) subject to the boundary conditions (11), (12) and (16) then leads 
to a characteristic equation for o in terms of the physical parameters of the problem. 
4. Solution and results. The derivation of the characteristic equation for o is not 
difficult but can be entirely avoided in the following manner. Observe that the dependence 
of the boundary condition (16) on 7, is more apparent than real; for, if we let 
° o1.oh” Otsy ¢ 


a = so that = +, 
Vv 71:0 a 


(17) 





then this third boundary condition can be written in the form 

a = ¢ {gl — 21[U"(1) — 2U(1)}} (18) 
and this relation no longer contains any physical parameters of the problem. It is, in 
fact, identical with the corresponding boundary condition for a self-gravitating sphere; 
the solution, in terms of a and g, must therefore be identical with the one given by Lamb 
and Chandrasekhar: 





a’. 2(1 — 1) | , q- 210. | ; 
_—_ = = { Q 
q* i ] q 1+ l + 1) q = 201 41/2(9) ’ (J ) 
where 
Qi4:2(Q = J 143/2(Q)/J 141/2(Q) (20) 


and the J’s are spherical Bessel functions. For g — ©, the limiting case of 
‘small viscosity’, we recover Lamb’s result 


1,5 = (Ll — 1)(21 + 1)vo/R? + to: . (21) 
Thus, the physical forces which maintain the spherical form enter only in the definition 
of the inviscid modes o;,. and these, as we have seen, can be completely eliminated by a 


suitable choice of dimensionless parameters. 
The characteristic equation (19) admits solutions of two types. For values of a’ 
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which exceed a certain critical value, damped oscillations will occur, but for values of a’ 
less than this critical value, two aperiodic modes of decay appear. The solution of Eq. 
(19) at this critical point for the principal mode J = 2 is [3] 


O2:ol° 7 = 3.69 and Ce-5/Cs:0 = 0.968. (22) 


For a drop of water surrounded by air (7, = 74 dynes/cm) this gives a radius R = 0.23 
mm. Drops larger than this critical radius will therefore tend to oscillate while smaller 
drops will be aperiodically damped. 
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ON THE DIFFRACTION OF AN ARBITRARY PULSE 
BY A WEDGE OR A CONE* 


BY LU TING (Polytechnic Institute of Brooklyn) 


Abstract. By virtue of Green’s Theorem, it is shown that for the diffraction of an 
arbitrary two-dimensional incident pulse by a wedge of angle yu, the ratio of the resultant 
velocity potential to the corresponding value of the incident pulse at the corner of 
the wedge at any instant is equal to 27/(2r — uw); and that for the diffraction of a three- 
dimensional pulse by a cone of solid angle w, the ratio at the vertex of the cone is equal 
to 4r/(4r — w). 

Two-dimensional space. The statement concerning diffraction of a pulse by a 
wedge is evidently true in the special case of an incident plane Heaviside pulse which 
was solved by Keller and Blank [1]. It therefore also follows for all incident pulses which 
are superpositions of plane Heaviside pulses, or limits of such superpositions. Since 
this includes all incident pulses it yields the preceding statement. However, these con- 
siderations depend upon knowing the exact solution in a special case which the follow- 
ing proof does not require.** 

Let ¢ = 0 be the instant at which the incident pulse ¢“” hits the corner of the wedge, 
which is located at the origin (x, = 0, x. = 0). Let h(x, , x.) and k(x, , x2) denote, re- 
spectively, ¢'’’ and ¢{" at an instant = — ¢, < 0 if the corner is absent. If G represents 
the domain in the x, — 2x, plane outside which both h and k vanish, then the origin must 
lie outside G. When the wedge is present, the region G lies outside the wedge if the 
incident disturbance ¢“” has not hit either side of the wedge at t = — t, < 0. Then the 
resultant disturbance ¢ at any instant t, > — ¢, fulfills the wave equation and the same 
initial conditions as that of ¢“’, i.e., in the region exterior to the wedge 


g(—to , 71 ,%2) = A(a, , 27) and 7] —ty ,% , 2s) = k(x, , Xe) (1) 


*Received Feb. 10, 1959; revised manuscript received May 27, 1959. 
**This paragraph is based on a private communication from Prof. J. B. Keller, Institute of Mathe- 
matical Sciences, New York University. 
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In addition ¢ fulfills the boundary condition d¢/dn = 0 on the two sides of the wedge, 
where 0/dn means normal derivative. 

To express the resultant velocity potential at the origin in terms of the initial data, 
it is necessary to reexamine the derivation of Volterra’s formula in order to take care of 
the boundary condition. The derivation begins from Green’s Theorem, which states [2]: 


a & / be PI ‘4 pk 4 
1 ov ov 0g 0g dg . 
| (, — =< 7 = No - dA + | ahs = a aA = @, (2) 
Ja dl OX; OX» Js ol OX, OX 
where v(é, 2; , %) = [(4, — 1° — af — x3)"; and n, , n , m2 are the # —, x, — and x, 


components of the unit vector 7 normal to the surface A in the ¢, 2 , Z2 space. Without 














Fia. 1. Diffraction by a wedge 


losing any generality the speed of sound is set equal to unity. The surface A consists 


of (Fig. 1 


(i) A*, which is the part of the characteristic cone 1/v = 0 lying in between the 
planes ¢ = — ¢,,¢ = ¢, — 6 and the boundary planes. 
(ii) Z, which is the part of the plane ¢ = — 4, , lying inside A*. 


(iii) ZZ, which is the part of the plane ¢ = ¢, — 4, lying inside A*. 

(iv) B, which is the part of the boundary planes, lving inside A* and between planes 
t= — ft, andt = #, — 6. 

The finite part of both integrals over A* vanishes as usual [2]. On the surface of the 
wedge the terms inside the brackets of the integrands become the normal derivative of 
v and ¢, respectively. The former is shown to be zero by performing the differentiation, 
and the latter equals zero according to the boundary condition. Therefore, the integral 
over B vanishes. 

As 6 — 0, the finite part of the first integral over 77 approaches (27 — yu) oft, , 0, 9) 
and the second integral approaches zero [2]. 

As a result, Eq. (2) together with the initial conditions, yields: 
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wa Ov 
(2r — pg(t, , 0,0) = | E —h- Jay, (3) 
J at 


where J is the region of G which lies inside the cone A*. 
From Volterra’s formula, the left side of the equation is equal to 2me" (t, , 0, 0); 
it follows that 
2 
g(t, , 0,0) = —— 
= 2 


~ 


e(t, , 0, 0). (4) 


3y the method of images it is shown that the relationship expressed in Eq. (4) is still 
valid, if the incident disturbance hits either or both sides of the wedge before it hits 
the corner [3]. 
Three-dimensional space. Let ¢ = 0 be the instant the incident pulse yg‘ hits the 
cone* whose vertex is located at the origin of the x, , z. — and x; — axes. In the absence 
of the cone, the incident potential at the origin can be expressed by Kirkchhoff’s formula 


(4] 
wenn: 1 [ ae” 1 arf ae ll ,, 
(t,0,0,0) = -— | Fle -(t) 2 ©. —1ar[ae" 1s 5 
. ir Js {\ Jon r | an | rant at Jj ei (5) 


where r is the distance from the origin, 
iis + mer? (i) 
d/dn means differentiation along the outward normal to the surface S, [y"’] = 


*Here the instant ¢ = 0 is defined differently from the previous case. For, even if g“ hit a portion 
of the surface of the cone excluding the vertex, the resultant potential in general cannot be obtained 
from yg“ by the method of images, as in the two-dimensional case. 


SURFACE 3S 





Fic. 2. Diffraction by a cone 








92 LU TING (Vol. XVIII, No. 1 


gy’ (x, , 2,23, — r) and [¢] is called the retarded value of ¢, and S denotes a surface 
whose shortest distance from the surface of the cone is greater than T > 0 (lig. 2). The 
reason for doing so will become evident later. 


To obtain the resultant velocity potential at the vertex of the cone, Green’s Theorem 


is applied: 





, ( P \ 
j fa] l 1 dl¢ 1 > , 
| $ le] — (4) — al, 1S + [ V [el] dV = 0, (6) 
J Si.te+lel on \Pr r on i 
where V is the volume bounded by the surface S, the two concentric spheres r = ¢ and 


r = R and the surface of the cone, and S, , o, = and B denote, respectively, the surfaces 


of the volume V. 


A repetition of the passages in the derivation of Kirkchhoff’s formula [4] gives 
; 1} ae]|\or , 1] a a 
| _{(@ ee “ —+-|=—|7dS = 0. (7) 
Js | rLdat on r Lon 
On the surface of the cone dr/dn = 0, while the boundary condition gives [dg/dn] = 0; 
therefore, the integral over B vanishes. 
Since ¢ = ¢'"’ for? < 0, the value of [¢] and its derivative on the pee r= Ris 
equal to that of [y**’] for R > ¢. For a finite value of ¢, the integral over = as R — = is 
identical with the corresponding integral, with ¢ replaced by ¢“. This iene tends 


to zero. This condition has been imposed on ¢"' 
As e — 0, the integral over ¢ approaches — (47 — w) ¢(/, 0, 0, 0) and Eq. (7) becomes 


i 


f 
-~ a oe: Io ar T de |) 
g(t, 0, 0, 0 <i | 4) < (4) a ep. 22a as. (8) 
lr —oaoJds\ mr r Lon ron | dt_]) 


Since the minimum distance between the surface S and the cone is greater than 7’, the 
value of [¢] and its derivatives on S will be identical with that of [g*'’] for any instant 


< 2T. From Eqs. (8) and (5), it yields 


in the derivation of Eq. (5). 


Ar : : 
g(t,0,0,0) = ———¢ “(t,0,0,0) for t < 27. Q) 


' 4r — w 


Since T is arbitrary, therefore, Eq. (9) is valid for any instant. 
This relationship can be interpreted intuitively as if the space were contracted by the 
)/4xr in the vicinity of the vertex so that the velocity potential 


cone by a factor of (47 — 
An analogous 


at the vertex is intensif fied by the reciprocal of the factor of contraction. 

interpretation ean be given to the two-dimensional problem. 
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INCORPORATION OF ERROR FUNCTION ABSORPTION 
COEFFICIENT IN TRANSPORT EQUATION* 


sy JIN H. CHIN** anp STUART W. CHURCHILL (The University of Michigan) 


1. Introduction. ‘The absorption of monochromatic radiation can be represented 
by the Bouguer-Beer law 


T=¢* (1) 


where 7 is the transmittance, x the path length and o the spectral absorption coefficient. 
However, absorption varies very rapidly in the region of a band and the average trans- 
mittance over a finite interval of wave length in such a region does not follow Eq. (1). 
The coefficient « cannot then be determined with an instrument of finite resolution. 
Various models have been proposed for an absorption band in order to derive expressions 
for the average transmittance. The Elsasser model [1] idealizes an absorption band as 
an infinite array of equally intense, equidistant lines, each of the Lorentz shape, and 
each with the same half width. For lines far apart relative to their half-width, this model 
vields an expression for the “average” transmittance at wave length \ which can be 
condensed to the following form 

T, = erfe (yx)'”. (2) 


The error function absorption coefficient y is related to the true spectral coefficient o as 
follows: 
\+Ad/2 
\1/2 —or 
erfe (yz) “" = = [ e** dy. (3) 
AX Jy-ans2 
The magnitude of A\ must be such that there is no contribution to the transmittance 
at the center of the interval from lines outside the interval. 
Values of y or the equivalent have been determined experimentally for water vapor, 
carbon dioxide and other gases. In contrast to o, y is a well-behaved function of wave 
length. For mixtures of two or more independently absorbing gases 


erfe (yx)'”? = erfe (y,x)'”- erfe (yox)'”- erfe (ysx)'? +++ . (4) 


2. Utilization of error function coefficient. The objective of this paper is to demon- 
strate how the error function absorption coefficient can be incorporated in a solution 
of the transport equation which describes radiative transfer. Solutions of the transport 
equation or its equivalent can be expressed in the following general form: 

i(R, Q) = Flor, 0,2 ---), (5) 
where 7 is the radiant intensity at any position denoted by the vector R in any direction 
denoted by the vector Q. The solution is assumed to be a function of oz, ¢,2 and other 
parameters, where x is some distance characterizing the position R and g¢, is the attenua- 
tion coefficient for single scattering. 

*Received May 13, 1959. 

**Now with the General Electric Co., Cincinnati, Ohio. 
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The parameters and variables other than o are assumed to be constant over a small 
interval AX. The “average”’ ornmeg over the interval AX is then 


A+ AX/ nh+Ad/2 


1 ae 
ax ae ., WR, @) dd = = | __, Plex) aa. (6) 


This \ — integration involving ox can be replaced by the following integral involving yz. 


ia rt” F(yxn) dn = 
xf F(cx) dX = == | ne (7) 
BA denen r J; “me ~ E") 

This substitution converts any solution of the transport equation or its equivalent 
for the monochromatic intensity in terms of ox and other variables and parameters into 
a new expression for the “average”? monochromatic intensity in terms of yx and the 
other variables and parameters, i.e., Eq. (5) is replaced by 


i(R, Q) = Ayr, o,x--: .), (8) 
where H(yx, o,2 -+-) is the function obtained from the integration on the right side 
of Eq. (7). This integration can be carried out analytically or graphically for a series of 
values of \ depending on the complexity of F(ox, ---). 

Alternative forms for the integration in terms of the parameter y are 
] ‘i . 1 fey u) du 2 ¢° Flyxr” + 1)] dv 
— F(ox) dz = if ,= | oes. (9 
Ad J) ails —u)” wT Jo v + | } 


The theorem represented by Eq. (7) and the two corollaries represented by Eq. (9) 
are established in the following section. 

3. Proof of the theorem and corollaries. Any function sectionally continuous in the 
interval (0, ~) may be expanded in terms of the normalized Laguerre polynomials [2] 


< 
f(é) a, a, (10) 


where L,(é) is the Laguerre polynomial of order »: 
, = f ’ , r ¥ °( ek ‘ ’ . ’ \ 
L,@) = & — ie - (te ) = (-—1) ie = Hp — (0 ] (11) 


and 


=f fee EO) a, (12) 


yp! 
On physical grounds F (ex) is a continuous or sectionally continuous function of ox 


and, therefore, can be expanded as follows: 


~—Gz 


= 4,(2ox)e 
F(ox) = F,(2ex) = >» a (13) 


Vv: 


v=0 


where F(2cz) is a function of 2ez obtained from the original function F(¢z). 
Therefore 


1 nr+Ad/2 co a 1 h+AX/2 
LP pas) ay = SELL Pe pean 
AD ad ox) dy > >! ‘a . L,(2ox)e an}, (14) 
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with 


—t 
-2f F@)L co dé. (15) 


at £7 


can be formed 


Since L,(2ex) is a polynomial of argument 2ox and degree v, L,(20x) e 
by summation of the products of multiples of powers of x with the differentials of e 


with respect to x. For instance, 


L,(2cx)e"" =e, (16) 
L,(2ex)e"* = (—2o0x + le 
d (17) 
=e” + (22) z “te * 
L.(2ox)e""" = [(20x)* — 4(2ox) + 2]e™”, 
= 4o’x’e* — 8exre * + 2c”, (18) 


2 d° or d —O2) 1 —oz 
= (42°) 2 (e**) + (82) y= (e ”*) + 2e”, 


and in general, 





) 


oz ‘ y - 
L,(2ox)¢ = (2x)” — (”*) + = (22)’ 





ov =) en a (Co) tee tole, (19) 
ax 


~ DT oY __ ay fe). 


(vy — n)! 'Pn! dx 


Therefore, 


1 A+ A4X/2 ; o ’ a,v! -_ F n L h+AX/2 ae 
Ar e ; F(cx) dX = p> > i> — aif (22) a-* Van [ e™ di? (20) 


The interchange of the order of integration and differentiation is allowable since 


ardt+Ar/2 


. a ec ax. 


Ad. 
exists by definition. 


Now erfe (yx)'”’ is the Laplace transform of f(¢) [3], where 


0 . Sate, 
f(t) =J 1/2 (21) 
\— at >, 


erfe (yx)' ‘= : [ Ys dt. (22) 
wT Jy — 
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Substituting ¢ = yy in the integral, yields 
co ~~ VEN 
erfe (yx)'”* = : od 
1 n(n — 1) 


From Eqs. (3) and (23) 
1 4+AX/2 1 2 a 7s 
an --/ ———173 dy. 
AX Jy-ans2 slits wJi n(n — 1)” ' 
From Eqs. (20) and (24), interchanging order of integration, differentiation, and sum- 


mation, yields 


pr+AX/2 


1 ee 7 sche es 
D hea - /F(ox) dy = } ee [v - = Pai 2 ») =~! n(n — 1)'”* in} (25) 





« ¥ v—n 


1 a,v! — . 
= a mE ————_,— (27)"™ yrn ‘ 
, rn(n — 1)” » 2, [(v — n)!fn! (22) dx’ (eo) dn 


dv 


But, it is evident from Eqs. (13) and (19) that 
~< a,y! wie 

———,— (2r)’~" —— (e**") = F(ynz). 

2 2 iw = n)tPat 22) gaan O°) = Flat) 


Therefore, 


F(ex) dX = 
Jx-An/2 on tJ, n(n — 1) 


\172 » 


ah+AX/2 a 
7 \ / 1 F(yzn) dy (27) 


and the theorem is proved. 

The proof of the corollaries is evident by substituting 7 = 1/u and n = v’ + 1in 
Eq. (27). 

4. Discussion. The incorporation of the error function coefficient in the formal 
solution of the transport equation immediately permits use of this relatively well- 
behaved coefficient in all existing solutions for light scattering, radiant transport, etc., 
subject only to the validity of the Elsasser model. The resulting ‘‘average’’ intensity, 
i, , can readily be integrated with respect to wave length for any given source distribution 
to find the total radiation over any finite interval of wave length. These procedures are 
illustrated in detail for a particular application in Ref. [4]. Attention is being given to 
equivalent theorems for other models for band absorption. 
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